hh_psc_alpha#

class brainpy.state.hh_psc_alpha(in_size, E_L=Quantity(-54.402, 'mV'), C_m=Quantity(100., 'pF'), g_Na=Quantity(12000., 'nS'), g_K=Quantity(3600., 'nS'), g_L=Quantity(30., 'nS'), E_Na=Quantity(50., 'mV'), E_K=Quantity(-77., 'mV'), t_ref=Quantity(2., 'ms'), tau_syn_ex=Quantity(0.2, 'ms'), tau_syn_in=Quantity(2., 'ms'), I_e=Quantity(0., 'pA'), V_m_init=Quantity(-65., 'mV'), Act_m_init=None, Inact_h_init=None, Act_n_init=None, spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', gsl_error_tol=0.001, name=None)#

NEST-compatible Hodgkin-Huxley neuron with alpha-shaped postsynaptic currents.

Current-based spiking neuron using the Hodgkin-Huxley formalism with voltage-gated sodium and potassium channels, leak conductance, alpha-function postsynaptic currents, threshold-and-local-maximum spike detection, and an explicit refractory period that suppresses spike emission only (subthreshold dynamics continue freely). Follows NEST models/hh_psc_alpha.{h,cpp} implementation with adaptive Runge-Kutta integration (RK45).

1. Mathematical Model

Membrane and ionic current dynamics:

The membrane potential evolves as

\[C_m \frac{dV_m}{dt} = -(I_{Na} + I_K + I_L) + I_{stim} + I_e + I_{syn,ex} + I_{syn,in}\]

where

\[\begin{split}I_{Na} &= g_{Na}\, m^3\, h\, (V_m - E_{Na}) \\ I_K &= g_K\, n^4\, (V_m - E_K) \\ I_L &= g_L\, (V_m - E_L)\end{split}\]

Gating variables \(m\) (Na activation), \(h\) (Na inactivation), \(n\) (K activation) obey first-order kinetics:

\[\frac{dx}{dt} = \alpha_x(V)(1 - x) - \beta_x(V)\,x\]

with voltage-dependent rate functions (voltage \(V\) in mV, rates in 1/ms):

\[\begin{split}\alpha_n &= \frac{0.01\,(V + 55)}{1 - e^{-(V+55)/10}}, \quad \beta_n = 0.125\,e^{-(V+65)/80} \\ \alpha_m &= \frac{0.1\,(V + 40)}{1 - e^{-(V+40)/10}}, \quad \beta_m = 4\,e^{-(V+65)/18} \\ \alpha_h &= 0.07\,e^{-(V+65)/20}, \quad \beta_h = \frac{1}{1 + e^{-(V+35)/10}}\end{split}\]

Alpha-function synaptic currents:

Each synapse type (excitatory/inhibitory) is modelled as a second-order linear system producing an alpha-shaped postsynaptic current:

\[\begin{split}\frac{dI_{syn}}{dt} &= dI_{syn} - \frac{I_{syn}}{\tau_{syn}} \\ \frac{d(dI_{syn})}{dt} &= -\frac{dI_{syn}}{\tau_{syn}}\end{split}\]

A spike arriving with weight \(w\) (in pA) adds \(w \cdot e / \tau_{syn}\) to \(dI_{syn}\), normalizing the peak current to \(w\) for \(w = 1\). Incoming spike weights are split by sign: positive weights drive excitatory state (\(dI_{syn,ex}\)), negative weights drive inhibitory state (\(dI_{syn,in}\)).

2. Spike Detection and Refractory Handling

A spike is detected when the membrane potential crosses 0 mV from below and a local maximum is detected (i.e., the potential starts decreasing). Formally, a spike is emitted when:

  1. refractory_step_count == 0 (not in refractory period), and

  2. V_m >= 0 mV (threshold crossing), and

  3. V_old > V_m (local maximum — potential is now falling).

Unlike integrate-and-fire models, no voltage reset occurs. The potassium current naturally repolarizes the membrane after a spike. During the refractory period \(t_{ref}\), spike emission is suppressed but all state variables continue evolving according to their differential equations.

3. Update Order Per Simulation Step

The update follows NEST’s exact order:

  1. Record pre-integration membrane potential (V_old).

  2. Integrate the full 8-dimensional ODE system \((V_m, m, h, n, dI_{ex}, I_{ex}, dI_{in}, I_{in})\) over one time step \([t, t+dt]\) using adaptive RK45 (Dormand-Prince).

  3. Add arriving synaptic spike inputs to \(dI_{syn,ex}\) and \(dI_{syn,in}\).

  4. Check spike condition: V_m >= 0 and V_old > V_m and r == 0.

  5. Update refractory counter and record spike time.

  6. Store buffered external stimulation current for the next step.

4. Numerical Integration

Uses AdaptiveRungeKuttaStep with method 'RKF45' to match NEST’s GSL RKF45 adaptive integrator. Default tolerance is gsl_error_tol=1e-3. All neurons are integrated simultaneously in a vectorized fashion.

5. Assumptions, Constraints, and Computational Implications

  • C_m > 0, g_Na >= 0, g_K >= 0, g_L >= 0, tau_syn_ex > 0, tau_syn_in > 0, and t_ref >= 0 are enforced at construction.

  • External current update(x=...) is buffered for one step, matching NEST ring-buffer semantics.

  • The adaptive RK45 integrator performs vectorized integration across all neurons simultaneously using JAX operations.

  • Spike detection uses a local maximum criterion rather than threshold crossing alone, matching biological action potential dynamics.

Parameters:
  • in_size (Size) – Population shape specification. All per-neuron parameters are broadcast to self.varshape derived from in_size.

  • E_L (ArrayLike, optional) – Leak reversal potential \(E_L\) in mV; scalar or array broadcastable to self.varshape. Determines resting potential. Default is -54.402 * u.mV.

  • C_m (ArrayLike, optional) – Membrane capacitance \(C_m\) in pF; broadcastable to self.varshape and strictly positive. Default is 100. * u.pF.

  • g_Na (ArrayLike, optional) – Sodium peak conductance \(g_{Na}\) in nS; broadcastable to self.varshape and non-negative. Default is 12000. * u.nS.

  • g_K (ArrayLike, optional) – Potassium peak conductance \(g_K\) in nS; broadcastable to self.varshape and non-negative. Default is 3600. * u.nS.

  • g_L (ArrayLike, optional) – Leak conductance \(g_L\) in nS; broadcastable to self.varshape and non-negative. Default is 30. * u.nS.

  • E_Na (ArrayLike, optional) – Sodium reversal potential \(E_{Na}\) in mV; broadcastable to self.varshape. Default is 50. * u.mV.

  • E_K (ArrayLike, optional) – Potassium reversal potential \(E_K\) in mV; broadcastable to self.varshape. Default is -77. * u.mV.

  • t_ref (ArrayLike, optional) – Absolute refractory period \(t_{ref}\) in ms; broadcastable to self.varshape and non-negative. Converted to integer step counts by ceil(t_ref / dt). Default is 2. * u.ms.

  • tau_syn_ex (ArrayLike, optional) – Excitatory alpha time constant \(\tau_{syn,ex}\) in ms; broadcastable to self.varshape and strictly positive. Default is 0.2 * u.ms.

  • tau_syn_in (ArrayLike, optional) – Inhibitory alpha time constant \(\tau_{syn,in}\) in ms; broadcastable to self.varshape and strictly positive. Default is 2. * u.ms.

  • I_e (ArrayLike, optional) – Constant injected current \(I_e\) in pA; scalar or array broadcastable to self.varshape. Default is 0. * u.pA.

  • V_m_init (ArrayLike, optional) – Initial membrane potential in mV; broadcastable to self.varshape. Default is -65. * u.mV.

  • Act_m_init (ArrayLike or None, optional) – Initial Na activation gating variable (dimensionless, range [0,1]); broadcastable to self.varshape. If None, initialized to equilibrium value at V_m_init. Default is None.

  • Inact_h_init (ArrayLike or None, optional) – Initial Na inactivation gating variable (dimensionless, range [0,1]); broadcastable to self.varshape. If None, initialized to equilibrium value at V_m_init. Default is None.

  • Act_n_init (ArrayLike or None, optional) – Initial K activation gating variable (dimensionless, range [0,1]); broadcastable to self.varshape. If None, initialized to equilibrium value at V_m_init. Default is None.

  • spk_fun (Callable, optional) – Surrogate spike nonlinearity used by get_spike() for differentiable spike generation. Default is braintools.surrogate.ReluGrad().

  • spk_reset (str, optional) – Reset policy inherited from Neuron. 'hard' applies stop-gradient to match NEST hard reset semantics. Default is 'hard'.

  • gsl_error_tol (float, optional) – Unitless local RKF45 error tolerance, strictly positive. Default is 1e-3.

  • name (str or None, optional) – Optional node name for debugging and visualization.

Parameter Mapping

Table 16 Parameter mapping to model symbols#

Parameter

Type / shape / unit

Default

Math symbol

Semantics

in_size

Size; scalar/tuple

required

Defines neuron population shape self.varshape.

E_L

ArrayLike, broadcastable to self.varshape (mV)

-54.402 * u.mV

\(E_L\)

Leak reversal potential (resting potential).

C_m

ArrayLike, broadcastable (pF), > 0

100. * u.pF

\(C_m\)

Membrane capacitance.

g_Na

ArrayLike, broadcastable (nS), >= 0

12000. * u.nS

\(g_{Na}\)

Sodium peak conductance.

g_K

ArrayLike, broadcastable (nS), >= 0

3600. * u.nS

\(g_K\)

Potassium peak conductance.

g_L

ArrayLike, broadcastable (nS), >= 0

30. * u.nS

\(g_L\)

Leak conductance.

E_Na

ArrayLike, broadcastable (mV)

50. * u.mV

\(E_{Na}\)

Sodium reversal potential.

E_K

ArrayLike, broadcastable (mV)

-77. * u.mV

\(E_K\)

Potassium reversal potential.

t_ref

ArrayLike, broadcastable (ms), >= 0

2. * u.ms

\(t_{ref}\)

Absolute refractory period duration.

tau_syn_ex

ArrayLike, broadcastable (ms), > 0

0.2 * u.ms

\(\tau_{syn,ex}\)

Excitatory alpha-kernel time constant.

tau_syn_in

ArrayLike, broadcastable (ms), > 0

2. * u.ms

\(\tau_{syn,in}\)

Inhibitory alpha-kernel time constant.

I_e

ArrayLike, broadcastable (pA)

0. * u.pA

\(I_e\)

Constant external current added every step.

V_m_init

ArrayLike, broadcastable (mV)

-65. * u.mV

Initial membrane potential.

Act_m_init

ArrayLike or None, dimensionless

None

Initial Na activation; None uses equilibrium at V_m_init.

Inact_h_init

ArrayLike or None, dimensionless

None

Initial Na inactivation; None uses equilibrium at V_m_init.

Act_n_init

ArrayLike or None, dimensionless

None

Initial K activation; None uses equilibrium at V_m_init.

spk_fun

Callable

ReluGrad()

Surrogate gradient function for spike generation.

spk_reset

str

'hard'

Reset mode; 'hard' stops gradient through reset.

gsl_error_tol

float, > 0

1e-3

Local absolute tolerance for the embedded RKF45 error estimate.

V#

Membrane potential \(V_m\). Shape: (*in_size,). Units: mV.

Type:

brainstate.HiddenState

m#

Na activation gating variable (dimensionless). Shape: (*in_size,). Range: [0, 1].

Type:

brainstate.HiddenState

h#

Na inactivation gating variable (dimensionless). Shape: (*in_size,). Range: [0, 1].

Type:

brainstate.HiddenState

n#

K activation gating variable (dimensionless). Shape: (*in_size,). Range: [0, 1].

Type:

brainstate.HiddenState

I_syn_ex#

Excitatory postsynaptic current \(I_{syn,ex}\). Shape: (*in_size,). Units: pA.

Type:

brainstate.ShortTermState

I_syn_in#

Inhibitory postsynaptic current \(I_{syn,in}\). Shape: (*in_size,). Units: pA.

Type:

brainstate.ShortTermState

dI_syn_ex#

Excitatory alpha-kernel derivative state. Shape: (*in_size,). Units: pA/ms.

Type:

brainstate.ShortTermState

dI_syn_in#

Inhibitory alpha-kernel derivative state. Shape: (*in_size,). Units: pA/ms.

Type:

brainstate.ShortTermState

I_stim#

One-step delayed external current buffer. Shape: (*in_size,). Units: pA.

Type:

brainstate.ShortTermState

refractory_step_count#

Remaining refractory steps. Shape: (*in_size,). Dtype: int32.

Type:

brainstate.ShortTermState

integration_step#

Persistent RKF45 substep size estimate (ms).

Type:

brainstate.ShortTermState

last_spike_time#

Time of most recent spike emission. Shape: (*in_size,). Units: ms. Updated to t + dt on spike emission.

Type:

brainstate.ShortTermState

Raises:

ValueError – If any of the following conditions are violated: - C_m <= 0 - t_ref < 0 - tau_syn_ex <= 0 or tau_syn_in <= 0 - g_Na < 0, g_K < 0, or g_L < 0

Notes

  • Unlike IAF models, the HH model does not reset the membrane potential after a spike. Repolarization occurs naturally through the potassium current.

  • During the refractory period, the neuron’s subthreshold dynamics continue to evolve freely; only spike emission is suppressed.

  • Spike weights are interpreted as current amplitudes (pA). Positive weights are excitatory; negative weights are inhibitory.

  • The adaptive RK45 integrator evaluates the ODE right-hand side multiple times per step, so computation cost scales with desired accuracy (controlled by gsl_error_tol).

  • Spike detection combines threshold crossing (0 mV) and local maximum detection, matching the biological action potential waveform.

References

See also

iaf_psc_alpha

Leaky integrate-and-fire with alpha-shaped PSCs.

hh_psc_alpha_clopath

HH neuron with Clopath voltage-based STDP.

hh_psc_alpha_gap

HH neuron with gap junction support.

Examples

Create a single Hodgkin-Huxley neuron and observe spiking behavior under constant current injection:

>>> import brainstate as bst
>>> import saiunit as u
>>> import brainpy.state as bps
>>> import matplotlib.pyplot as plt
>>> # Initialize simulation context
>>> bst.environ.set(dt=0.1 * u.ms)
>>> # Create neuron
>>> neuron = bps.hh_psc_alpha(in_size=1, I_e=500. * u.pA)
>>> neuron.init_all_states()
>>> # Run simulation
>>> times = []
>>> voltages = []
>>> for _ in range(2000):  # 200 ms
...     neuron.update()
...     times.append(float(bst.environ.get('t') / u.ms))
...     voltages.append(float(neuron.V.value / u.mV))
>>> # Plot results
>>> plt.plot(times, voltages)
>>> plt.xlabel('Time (ms)')
>>> plt.ylabel('Membrane potential (mV)')
>>> plt.title('Hodgkin-Huxley neuron spiking')
>>> plt.show()
get_spike(V=None)[source]#

Generate differentiable spike output using surrogate gradient.

Applies the surrogate spike function to the membrane potential scaled relative to the 0 mV threshold. This enables gradient-based learning through the spike generation process.

Parameters:

V (ArrayLike or None, optional) – Membrane potential in mV. If None, uses self.V.value. Shape must broadcast with self.varshape.

Returns:

Differentiable spike signal with shape (*in_size,). Typically near 0 for subthreshold, near 1 for suprathreshold.

Return type:

ArrayLike

Notes

The spike threshold for HH neurons is 0 mV. The input voltage is scaled relative to this threshold before applying the surrogate function.

init_state(**kwargs)[source]#

Initialize all state variables.

Sets initial values for membrane potential, gating variables, synaptic currents, refractory counters, and buffers. If gating variable initial values are not explicitly provided, they are computed at equilibrium for the given initial membrane potential.

Parameters:

**kwargs (dict) – Additional keyword arguments (unused, for compatibility).

Notes

State variables initialized:
  • V: membrane potential (from V_m_init)

  • m, h, n: gating variables (from Act_m_init, Inact_h_init, Act_n_init if provided; otherwise computed at equilibrium for V_m_init)

  • I_syn_ex, I_syn_in, dI_syn_ex, dI_syn_in: synaptic states (initialized to zero)

  • I_stim: external current buffer (initialized to zero)

  • refractory_step_count: refractory countdown (initialized to zero)

  • integration_step: persistent RKF45 substep size

  • last_spike_time: spike time record (initialized to -1e7 ms)

update(x=Quantity(0., 'pA'))[source]#

Update neuron state for one simulation step.

Integrates the full Hodgkin-Huxley dynamics over one time step \(dt\), applies synaptic inputs, detects spikes using threshold-and-local-maximum criterion, updates refractory state, and buffers external current for the next step. Follows NEST hh_psc_alpha update order exactly.

Update Order:

  1. Record pre-integration membrane potential (V_old).

  2. Integrate the 8-dimensional ODE system \((V_m, m, h, n, dI_{ex}, I_{ex}, dI_{in}, I_{in})\) over \([t, t+dt]\) using adaptive RK45 (Dormand-Prince).

  3. Add arriving synaptic spike inputs to \(dI_{syn,ex}\) and \(dI_{syn,in}\).

  4. Check spike condition: refractory_step_count == 0 and V_m >= 0 and V_old > V_m.

  5. Update refractory counter and record spike time.

  6. Store buffered external stimulation current x for next step.

Parameters:

x (ArrayLike, optional) – External stimulation current input in pA (in addition to I_e). Shape must broadcast with (*in_size,). Default is 0. * u.pA.

Returns:

Differentiable spike output with shape (*in_size,). Generated by applying self.spk_fun to the spike condition. Near 1 when spike detected, near 0 otherwise.

Return type:

ArrayLike

Notes

  • The external current x is buffered for one step via I_stim, matching NEST’s ring-buffer semantics. Current provided at step \(n\) affects dynamics at step \(n+1\).

  • Spike weights are collected via sum_delta_inputs(0*pA) and split by sign: positive weights drive excitatory state, negative weights drive inhibitory state.

  • During the refractory period, all state variables evolve freely; only spike emission is suppressed.

  • Spike detection combines threshold crossing (0 mV) and local maximum detection (V_old > V_m) to match biological action potential characteristics.