Hodgkin-Huxley
The Hodgkin-Huxley neuron, whose sodium and potassium channels produce action potentials.
- Name:
hodgkin_huxley - Field:
biology/neuroscience - Formalism: ordinary differential equations
- References: Hodgkin and Huxley 1952; Rinzel and Miller 1980; Trzaska 2023; arXiv:2112.12147
Factory
@register_model(
field='biology/neuroscience',
tags=('neuron', 'excitability', 'spiking', 'conductance'),
references=('Hodgkin and Huxley 1952', 'Rinzel and Miller 1980', 'Trzaska 2023', 'arXiv:2112.12147'),
defaults={'until': 50.0, 'dt': 0.1},
)
def hodgkin_huxley(
v: float = -65.0,
*,
I: float = 10.0,
C: float = 1.0,
g_na: float = 120.0,
g_k: float = 36.0,
g_l: float = 0.3,
) -> System:
"""Return the Hodgkin-Huxley neuron, whose sodium and potassium channels produce action potentials.
The gates start at their steady state at the initial potential, and the current is switched on at the start. A
constant current above about 6.2 µA/cm² makes the neuron fire periodically, and a weaker one gives at most a few
spikes before it settles.
Args:
v:
The initial membrane potential in mV.
I:
The applied current in µA/cm².
C:
The membrane capacitance in µF/cm².
g_na:
The maximal sodium conductance in mS/cm².
g_k:
The maximal potassium conductance in mS/cm².
g_l:
The leak conductance in mS/cm².
Returns:
system:
The model, with the potential `neuron.v`, the gates `neuron.n`, `neuron.m`, and `neuron.h`, and its spikes
as `spike` events.
Examples:
>>> from skcomplex.models.biology.neuroscience import hodgkin_huxley
>>> round(hodgkin_huxley().elements[0].h, 4)
0.5961
"""
gates = _resting_gates(v)
return System(
elements=[
HodgkinHuxleyElement(
'neuron',
v=v,
I=I,
C=C,
g_na=g_na,
g_k=g_k,
g_l=g_l,
n=gates['n'],
m=gates['m'],
h=gates['h'],
),
],
interactions=[
RateInteraction('applied_current', target='neuron', state='v', rate=lambda s: s.target.I, divisor='C'),
RateInteraction(
'sodium_current',
target='neuron',
state='v',
rate=lambda s: -s.target.g_na * s.target.m**3 * s.target.h * (s.target.v - s.target.e_na),
divisor='C',
),
RateInteraction(
'potassium_current',
target='neuron',
state='v',
rate=lambda s: -s.target.g_k * s.target.n**4 * (s.target.v - s.target.e_k),
divisor='C',
),
RateInteraction(
'leak_current',
target='neuron',
state='v',
rate=lambda s: -s.target.g_l * (s.target.v - s.target.e_l),
divisor='C',
),
*[
RateInteraction(f'{gate}_gating', target='neuron', state=gate, rate=_gate_kinetics(gate))
for gate in ('n', 'm', 'h')
],
CrossingInteraction('spike', target='neuron', guard=lambda s: s.target.v, direction=1),
],
)
Components
class HodgkinHuxleyElement(Element):
"""Squid giant axon membrane of Hodgkin and Huxley, with sodium, potassium, and leak currents.
The potential follows `C v' = I - g_Na m³ h (v - E_Na) - g_K n⁴ (v - E_K) - g_L (v - E_L)`, and each gate `x` of
`n`, `m`, and `h` relaxes as `x' = alpha_x(v) (1 - x) - beta_x(v) x`, in the modern convention with rest near -65
mV. Upward crossings of 0 mV are recorded as `spike` events. Units are mV, ms, mS/cm², and µA/cm². The element holds
the states and parameters and has no laws of its own: the current, gating, and crossing interactions of
`hodgkin_huxley` read its parameters.
Args:
label:
The label of the neuron.
visual:
How the element is drawn, or `None` for `default_visual_`.
v:
The initial membrane potential.
n:
The initial potassium activation.
m:
The initial sodium activation.
h:
The initial sodium inactivation.
C:
The membrane capacitance.
g_na:
The maximal sodium conductance.
g_k:
The maximal potassium conductance.
g_l:
The leak conductance.
e_na:
The sodium reversal potential.
e_k:
The potassium reversal potential.
e_l:
The leak reversal potential.
I:
The applied current.
Examples:
>>> from skcomplex.models.biology.neuroscience import HodgkinHuxleyElement
>>> HodgkinHuxleyElement('axon', I=10.0).g_na
120.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
v: Annotated[float, StateSpec()] = -65.0,
n: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.3177,
m: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0529,
h: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5961,
C: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
g_na: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 120.0,
g_k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 36.0,
g_l: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
e_na: Annotated[float, ParameterSpec()] = 50.0,
e_k: Annotated[float, ParameterSpec()] = -77.0,
e_l: Annotated[float, ParameterSpec()] = -54.387,
I: Annotated[float, ParameterSpec()] = 10.0,
) -> None:
super().__init__(label, visual=visual)
self.v = v
self.n = n
self.m = m
self.h = h
self.C = C
self.g_na = g_na
self.g_k = g_k
self.g_l = g_l
self.e_na = e_na
self.e_k = e_k
self.e_l = e_l
self.I = I
class RateInteraction(Interaction):
"""Rate `rate(ns)` and noise `noise(ns) dW` added to the highest time derivative of a state of a target.
On a state of order 1 the rate is its first derivative, and on a state of order `k` it adds to the `k`-th
derivative, so rates on one state superpose. With a divisor `d` the rate and the noise are divided by `d`, the name
of a state or parameter of the target, such as an inertia, a friction, or a capacitance, which backends can read.
A rate, a noise, or a divisor function takes the namespace of the interaction, with the namespaces `target` and
`sources` of the bound elements, the time `t`, and the symbols of the parameters a subclass declares. A number
stands for a constant, and a number on a vector state applies to every component. The noise is a Wiener process per
component and member by default, a process of a `WienerProcessSpec` to share or correlate it across states and
members, or a fractional Brownian motion of a Hurst index.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state the rate changes.
state:
The name of the state of the target that the rate changes.
rate:
A function of the namespace that returns the rate, or a number.
sources:
The paths of the elements or groups the rate and the noise read but do not change.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
divisor:
The name of the state or parameter of the target that divides the rate and the noise, a function of the
namespace of the target that returns the divisor, such as a moment of inertia `m l^2`, or `None`.
process:
The Wiener processes that drive the noise, or `None` for one independent process per component and member.
index:
The process of `process` that drives the noise, or `None` for every process, one per component of a vector
state.
convention:
The stochastic calculus of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.
hurst:
The Hurst index in `(0, 1)` of a fractional Brownian motion that drives the noise, a number or a function
of the namespace that returns a parameter, or `None` for a Wiener process.
conserved:
Whether the noise is the divergence of a random flux on the links of the lattice of a field, so it keeps the
total of the state.
Raises:
ValidationError:
The state is not a name, a function is neither callable nor a number, the convention is unknown, the Hurst
index is outside `(0, 1)`, or the noise options cannot combine.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RateInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> decay = RateInteraction('decay', target='stock', state='n', rate=lambda s: -0.5 * s.target.n)
>>> system = System(elements=[CompartmentElement('stock', n=8.0)], interactions=[decay])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('stock.n')[-1]), 3)
2.943
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
rate: Callable[[SimpleNamespace], object] | float = 0.0,
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
noise: Callable[[SimpleNamespace], object] | float | None = None,
divisor: str | Callable[[SimpleNamespace], object] | None = None,
process: Annotated[WienerProcessSpec | None, RoleSpec(noise=True)] = None,
index: int | None = None,
convention: str = 'ito',
hurst: Callable[[SimpleNamespace], object] | float | None = None,
conserved: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.rate = rate
self.sources = sources
self.noise = noise
self.divisor = divisor
self.process = process
self.index = index
self.convention = convention
self.hurst = hurst
self.conserved = conserved
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} changes the state {self.state!r}, which is not a name of a '
'state.'
)
raise ValidationError(error_msg)
for name in ('rate', 'noise'):
value = getattr(self, name)
if not (value is None or callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
if not (
self.divisor is None
or callable(self.divisor)
or (isinstance(self.divisor, str) and self.divisor.isidentifier())
):
error_msg = (
f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, or a '
'function of the namespace of the target.'
)
raise ValidationError(error_msg)
def _check_noise(self) -> None:
"""Raise an error for an unknown convention, a Hurst index outside `(0, 1)`, or noise options that conflict."""
if self.convention not in CONVENTIONS:
error_msg = f'Unknown noise convention {self.convention!r}. Use one of {", ".join(CONVENTIONS)}.'
raise ValidationError(error_msg)
if not (
self.hurst is None or callable(self.hurst) or (isinstance(self.hurst, int | float) and 0 < self.hurst < 1)
):
error_msg = f'The Hurst index of {self.label!r} is {self.hurst!r}. Give a number in (0, 1) or a function.'
raise ValidationError(error_msg)
if self.index is not None and self.process is None:
error_msg = f'{self.label!r} picks the process {self.index!r} but binds no processes. Give `process`.'
raise ValidationError(error_msg)
options = [
name
for name, given in (('hurst', self.hurst), ('process', self.process), ('conserved', self.conserved))
if given
]
if len(options) > 1:
error_msg = f'The noise of {self.label!r} takes {" and ".join(options)}. Give at most one of them.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for a state, functions, or noise options the interaction cannot take.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is neither callable nor a number, the
convention is unknown, the Hurst index is outside `(0, 1)`, or the noise options cannot combine.
"""
super()._check()
self._check_functions()
self._check_noise()
def _evaluate(self, function: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
"""Return the value of a function of the namespace, or a number as it is."""
return function(namespace) if callable(function) else function
def _broadcast(self, value: object, state: object) -> object:
"""Return a value with one component per component of a vector state, repeating a number."""
if isinstance(state, FreeModuleElement) and not isinstance(value, FreeModuleElement | list | tuple):
return vector(SR, [value] * len(state))
if isinstance(value, list | tuple):
return vector(SR, value)
return value
def _divide(self, value: object, namespace: SimpleNamespace) -> object:
"""Return a value divided by the divisor of the target, or the value without a divisor."""
if self.divisor is None:
return value
divisor = (
self.divisor(self._get_member(namespace, 'target'))
if callable(self.divisor)
else getattr(self._get_member(namespace, 'target'), self.divisor)
)
return operator.truediv(value, divisor)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the rate before the divisor, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.rate, namespace)
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the noise amplitude before the divisor, the hook a subclass overrides, or `None` for no noise."""
return None if self.noise is None else self._evaluate(self.noise, namespace)
def _get_noise_law(self, state: object, amplitude: object, namespace: SimpleNamespace) -> BasePart:
"""Return the noise law of an amplitude, by the process, the convention, and the Hurst index."""
if self.hurst is not None:
hurst = self._evaluate(self.hurst, namespace)
return FractionalNoiseLaw(state, amplitude, hurst=hurst, convention=self.convention)
process = None
if self.process is not None:
handle = self._get_member(namespace, 'process').dW
process = handle if self.index is None else handle[self.index]
return NoiseLaw(state, amplitude, convention=self.convention, process=process, conserved=self.conserved)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the state and its noise law, both divided by the divisor.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the target and the sources.
Returns:
laws:
The rate law, and the noise law when the interaction has noise.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
rate = self._divide(self._broadcast(self._get_rate(namespace), state), namespace)
laws: list[BasePart] = [RateLaw(state, rate)]
noise = self._get_noise(namespace)
if noise is not None:
laws.append(self._get_noise_law(state, self._divide(noise, namespace), namespace))
return laws
class CrossingInteraction(Interaction):
"""Event at the crossing of zero by a guard, which sets states of the target to new values or ends the run.
The guard and the new values are functions of the namespace, which holds `namespace.target`, the states and
parameters of the target, and the time `namespace.t`. A guard that reads the members of a group gives every member
its own event, which sets the states of that member alone. Solvers locate the crossing within their tolerance, and
the trajectory records every crossing among its events, under the label of the interaction and with the member.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states the guard reads and the event sets.
guard:
The function `guard(namespace)` of the expression whose crossing of zero triggers the event.
reset:
A mapping from the name of a state of the target to the function `f(namespace)` of its new value, applied
all at once, or `None` to keep the states.
direction:
`1` for crossings from negative to positive, `-1` for the opposite, and `0` for both.
terminal:
Whether the crossing ends the run.
Raises:
ValidationError:
The direction is not -1, 0, or 1, or the guard or a new value is not a function.
Examples:
>>> from typing import Annotated
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.interactions import CrossingInteraction
>>> from skcomplex.core.laws import RateLaw
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class ClockElement(Element):
... def __init__(self, label, x: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.x = x
...
... def _get_laws(self, s):
... return [RateLaw(s.x, 1.0)]
>>> wrap = CrossingInteraction(
... 'wrap', target='clock', guard=lambda s: s.target.x - 1, reset={'x': lambda s: 0}, direction=1
... )
>>> system = System(elements=[ClockElement('clock')], interactions=[wrap])
>>> trajectory = Simulation(system, dt=0.5).check().prepare().run(until=2.2).trajectory_
>>> [round(event.time, 6) for event in trajectory.events]
[1.0, 2.0]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
guard: Callable[[SimpleNamespace], object],
reset: Mapping[str, Callable[[SimpleNamespace], object]] | None = None,
direction: int = 0,
terminal: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.guard = guard
self.reset = reset
self.direction = direction
self.terminal = terminal
def _check(self) -> None:
"""Raise an error for a direction that is not -1, 0, or 1, or a guard or a new value that is not a function.
Raises:
ValidationError:
The direction is not -1, 0, or 1, or the guard or a new value is not a function.
"""
super()._check()
if self.direction not in DIRECTIONS:
error_msg = f'The direction of crossing {self.label!r} is {self.direction!r}. Use -1, 0, or 1.'
raise ValidationError(error_msg)
if not callable(self.guard) or not all(callable(function) for function in (self.reset or {}).values()):
error_msg = f'The guard and the new values of crossing {self.label!r} must be functions of the namespace.'
raise ValidationError(error_msg)
@staticmethod
def _get_law_value(value: object) -> object:
"""Return what a function of the namespace returned: a function as it is, or a symbolic expression."""
return value if BaseLaw._is_function(value) else SR(value)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the crossing law of the guard with the new values of the states.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target`.
Returns:
laws:
The crossing law.
"""
reset = {
getattr(self._get_member(namespace, 'target'), name): self._get_law_value(f(namespace))
for name, f in (self.reset or {}).items()
}
return [
CrossingLaw(
self._get_law_value(self.guard(namespace)),
reset or None,
self.direction,
self.terminal,
label=self.label,
),
]
Equations
neuron.v: diff(v_neuron(t), t) == (e_na_neuron - v_neuron(t))*g_na_neuron*h_neuron(t)*m_neuron(t)^3/C_neuron + (e_k_neuron - v_neuron(t))*g_k_neuron*n_neuron(t)^4/C_neuron + (e_l_neuron - v_neuron(t))*g_l_neuron/C_neuron + I_neuron/C_neuron
neuron.n: diff(n_neuron(t), t) == -0.125*e^(-1/80*v_neuron(t) - 13/16)*n_neuron(t) + (n_neuron(t) - 1)*(0.01*v_neuron(t) + 0.55)/(e^(-1/10*v_neuron(t) - 11/2) - 1)
neuron.m: diff(m_neuron(t), t) == -4*e^(-1/18*v_neuron(t) - 65/18)*m_neuron(t) + (m_neuron(t) - 1)*(0.1*v_neuron(t) + 4.0)/(e^(-1/10*v_neuron(t) - 4) - 1)
neuron.h: diff(h_neuron(t), t) == -0.07*(h_neuron(t) - 1)*e^(-1/20*v_neuron(t) - 13/4) - h_neuron(t)/(e^(-1/10*v_neuron(t) - 7/2) + 1)