hh_cond_exp_traub

class brainpy.state.hh_cond_exp_traub(in_size, E_L=Quantity(-60., 'mV'), C_m=Quantity(200., 'pF'), g_Na=Quantity(20000., 'nS'), g_K=Quantity(6000., 'nS'), g_L=Quantity(10., 'nS'), E_Na=Quantity(50., 'mV'), E_K=Quantity(-90., 'mV'), V_T=Quantity(-63., 'mV'), E_ex=Quantity(0., 'mV'), E_in=Quantity(-80., 'mV'), t_ref=Quantity(2., 'ms'), tau_syn_ex=Quantity(5., 'ms'), tau_syn_in=Quantity(10., 'ms'), I_e=Quantity(0., 'pA'), V_m_init=None, Act_m_init=None, Inact_h_init=None, Act_n_init=None, gsl_error_tol=0.001, spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', name=None)

NEST-compatible hh_cond_exp_traub neuron model.

Hodgkin-Huxley model for Brette et al. (2007) review, based on Traub and Miles (1991) hippocampal pyramidal cell model.

This is a modified Hodgkin-Huxley neuron model specifically developed for the Brette et al. (2007) simulator review, based on a model of hippocampal pyramidal cells by Traub and Miles (1991). Key differences from the original Traub-Miles model:

• This is a point neuron, not a compartmental model.

• Only I_Na and I_K ionic currents are included (no calcium dynamics), with simplified I_K dynamics giving three gating variables instead of eight.

• Incoming spikes induce an instantaneous conductance change followed by exponential decay (conductance-based synapses), not activation over time.

Parameters:

• in_size (int, tuple of int) – Population shape (number of neurons or spatial dimensions).

• E_L (ArrayLike, default -60 mV) – Leak reversal potential. Must be finite.

• C_m (ArrayLike, default 200 pF) – Membrane capacitance. Must be strictly positive.

• g_Na (ArrayLike, default 20000 nS) – Sodium peak conductance. Must be non-negative.

• g_K (ArrayLike, default 6000 nS) – Potassium peak conductance. Must be non-negative.

• g_L (ArrayLike, default 10 nS) – Leak conductance. Must be non-negative.

• E_Na (ArrayLike, default 50 mV) – Sodium reversal potential. Must be finite.

• E_K (ArrayLike, default -90 mV) – Potassium reversal potential. Must be finite.

• V_T (ArrayLike, default -63 mV) – Voltage offset for gating dynamics. Shifts the effective threshold to approximately V_T + 30 mV.

• E_ex (ArrayLike, default 0 mV) – Excitatory synaptic reversal potential. Must be finite.

• E_in (ArrayLike, default -80 mV) – Inhibitory synaptic reversal potential. Must be finite.

• t_ref (ArrayLike, default 2 ms) – Duration of refractory period. Must be non-negative. Traub and Miles used 3 ms; NEST default is 2 ms.

• tau_syn_ex (ArrayLike, default 5 ms) – Excitatory synaptic time constant. Must be strictly positive.

• tau_syn_in (ArrayLike, default 10 ms) – Inhibitory synaptic time constant. Must be strictly positive.

• I_e (ArrayLike, default 0 pA) – Constant external input current. Can be positive or negative.

• V_m_init (ArrayLike, optional) – Initial membrane potential. If None, defaults to E_L.

• Act_m_init (ArrayLike, optional) – Initial sodium activation gating variable (0 <= m <= 1). If None, computed from equilibrium at V_m_init.

• Inact_h_init (ArrayLike, optional) – Initial sodium inactivation gating variable (0 <= h <= 1). If None, computed from equilibrium at V_m_init.

• Act_n_init (ArrayLike, optional) – Initial potassium activation gating variable (0 <= n <= 1). If None, computed from equilibrium at V_m_init.

• gsl_error_tol (ArrayLike) – Unitless local RKF45 error tolerance, broadcastable and strictly positive.

• spk_fun (Callable, default braintools.surrogate.ReluGrad()) – Surrogate spike function for differentiable spike generation.

• spk_reset (str, default 'hard') – Reset mode (‘hard’ or ‘soft’). Note: HH models do not reset voltage after spikes; this parameter affects gradient computation only.

• name (str, optional) – Name of the neuron population.

V

Membrane potential with shape (*in_size,) in mV.

Type:

brainstate.HiddenState

m

Sodium activation gating variable (0 <= m <= 1), shape (*in_size,).

Type:

brainstate.HiddenState

h

Sodium inactivation gating variable (0 <= h <= 1), shape (*in_size,).

Type:

brainstate.HiddenState

n

Potassium activation gating variable (0 <= n <= 1), shape (*in_size,).

Type:

brainstate.HiddenState

g_ex

Excitatory synaptic conductance in nS, shape (*in_size,).

Type:

brainstate.HiddenState

g_in

Inhibitory synaptic conductance in nS, shape (*in_size,).

Type:

brainstate.HiddenState

I_stim

Stimulation current buffer in pA, shape (*in_size,).

Type:

brainstate.ShortTermState

refractory_step_count

Refractory countdown in grid 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 in ms, shape (*in_size,).

Type:

brainstate.ShortTermState

Raises:

ValueError – If C_m <= 0, t_ref < 0, tau_syn_ex <= 0, or tau_syn_in <= 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, subthreshold dynamics continue to evolve freely; only spike emission is suppressed.

• Synaptic spike weights are interpreted in conductance units (nS). Positive weights drive excitatory synapses; negative weights drive inhibitory synapses (sign is flipped, i.e. g_in += |w|).

• The numerical integration uses an adaptive RKF45 (Runge-Kutta-Fehlberg) integrator implemented in JAX with unit-aware arithmetic via saiunit. This is equivalent to NEST’s GSL RKF45 implementation for numerical correspondence.

Mathematical Formulation

1. Membrane and Ionic Current Dynamics

The membrane potential evolves as:

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

where the currents are:

\[\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) \\ I_{syn,ex} &= g_{ex}\, (V_m - E_{ex}) \\ I_{syn,in} &= g_{in}\, (V_m - E_{in})\end{split}\]

2. Channel Gating Variables

Gating variables \(m\), \(h\), \(n\) obey first-order kinetics:

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

with Traub-Miles rate functions using shifted voltage \(V = V_m - V_T\) (voltage in mV, rates in 1/ms):

\[\begin{split}\alpha_n &= \frac{0.032\,(15 - V)}{e^{(15 - V)/5} - 1}, \quad \beta_n = 0.5\,e^{(10 - V)/40} \\ \alpha_m &= \frac{0.32\,(13 - V)}{e^{(13 - V)/4} - 1}, \quad \beta_m = \frac{0.28\,(V - 40)}{e^{(V - 40)/5} - 1} \\ \alpha_h &= 0.128\,e^{(17 - V)/18}, \quad \beta_h = \frac{4}{1 + e^{(40 - V)/5}}\end{split}\]

The voltage offset \(V_T\) (default -63 mV) shifts the effective threshold to approximately -50 mV.

3. Exponential Conductance Synapses

Synaptic conductances decay exponentially:

\[\begin{split}\frac{dg_{ex}}{dt} &= -g_{ex} / \tau_{syn,ex} \\ \frac{dg_{in}}{dt} &= -g_{in} / \tau_{syn,in}\end{split}\]

A presynaptic spike with weight \(w\) causes an instantaneous conductance jump:

• \(w > 0\) – \(g_{ex} \leftarrow g_{ex} + w\)

• \(w < 0\) – \(g_{in} \leftarrow g_{in} + |w|\)

4. Spike Detection

A spike is emitted when all three conditions are satisfied:

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

2. V_m >= V_T + 30 mV (threshold crossing), and

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

Unlike integrate-and-fire models, no voltage reset occurs – the potassium current naturally repolarizes the membrane.

Warning

To avoid multiple spikes during the falling flank of a spike, it is essential to choose a sufficiently long refractory period. Traub and Miles used \(t_{ref} = 3\) ms, while the default here is \(t_{ref} = 2\) ms (matching NEST).

5. Numerical Integration

NEST uses GSL RKF45 (Runge-Kutta-Fehlberg 4/5) with adaptive step-size control. This implementation uses an adaptive RKF45 integrator implemented in JAX with unit-aware arithmetic via saiunit, matching NEST’s integration approach for numerical correspondence.

The ODE system is 6-dimensional per neuron: \([V_m, m, h, n, g_{ex}, g_{in}]\).

Parameter Mapping

The following table shows the correspondence between brainpy.state parameters and NEST/mathematical notation:

Parameter	Default	Math equivalent	Description
in_size	(required)	–	Population shape
E_L	-60 mV	\(E_L\)	Leak reversal potential
C_m	200 pF	\(C_m\)	Membrane capacitance
g_Na	20000 nS	\(g_{Na}\)	Sodium peak conductance
g_K	6000 nS	\(g_K\)	Potassium peak conductance
g_L	10 nS	\(g_L\)	Leak conductance
E_Na	50 mV	\(E_{Na}\)	Sodium reversal potential
E_K	-90 mV	\(E_K\)	Potassium reversal potential
V_T	-63 mV	\(V_T\)	Voltage offset for gating dynamics
E_ex	0 mV	\(E_{ex}\)	Excitatory synaptic reversal potential
E_in	-80 mV	\(E_{in}\)	Inhibitory synaptic reversal potential
t_ref	2 ms	\(t_{ref}\)	Duration of refractory period
tau_syn_ex	5 ms	\(\tau_{syn,ex}\)	Excitatory synaptic time constant
tau_syn_in	10 ms	\(\tau_{syn,in}\)	Inhibitory synaptic time constant
I_e	0 pA	\(I_e\)	Constant external input current
V_m_init	None	–	Initial V_m (None -> E_L)
Act_m_init	None	–	Initial Na activation (None -> equilibrium)
Inact_h_init	None	–	Initial Na inactivation (None -> equilibrium)
Act_n_init	None	–	Initial K activation (None -> equilibrium)
gsl_error_tol	1e-3	–	Local RKF45 error tolerance
spk_fun	ReluGrad()	–	Surrogate spike function
spk_reset	'hard'	–	Reset mode

Examples

>>> import brainstate as bst
>>> import saiunit as u
>>> from brainpy_state import hh_cond_exp_traub
>>>
>>> # Create a population of 100 Traub HH neurons
>>> neurons = hh_cond_exp_traub(100)
>>> neurons.init_all_states()
>>>
>>> # Run a simulation with constant current injection
>>> with bst.environ.context(dt=0.1*u.ms):
...     for i in range(1000):
...         spikes = neurons.update(I_e=200*u.pA)

>>> # Compare with NEST default parameters
>>> import nest
>>> nest_neuron = nest.Create('hh_cond_exp_traub')
>>> nest.GetStatus(nest_neuron, ['V_m', 'E_L', 'C_m', 'g_Na', 'g_K'])
[(-60.0, -60.0, 200.0, 20000.0, 6000.0)]
>>>
>>> # Match in brainpy.state
>>> bp_neuron = hh_cond_exp_traub(1, E_L=-60*u.mV, C_m=200*u.pF,
...                               g_Na=20000*u.nS, g_K=6000*u.nS)

References

[1]

Brette R et al. (2007). Simulation of networks of spiking neurons: A review of tools and strategies. Journal of Computational Neuroscience 23:349-98. DOI: https://doi.org/10.1007/s10827-007-0038-6

[2]

Traub RD and Miles R (1991). Neuronal networks of the hippocampus. Cambridge University Press, Cambridge UK.

[3]

ModelDB entry: http://modeldb.yale.edu/83319

See also

hh_psc_alpha

Hodgkin-Huxley with alpha-shaped postsynaptic currents.

iaf_cond_exp

Leaky integrate-and-fire with conductance-based synapses.

get_spike(V=None)

Compute differentiable spike output using surrogate gradient function.

Applies the surrogate spike function to the membrane potential. This is used for gradient-based learning; actual spike detection in the update method uses discrete threshold crossing logic (V >= V_T + 30 and local maximum).

Parameters:

V (ArrayLike, optional) – Membrane potential in mV, shape (*in_size,) or (batch_size, *in_size). If None, uses the current state self.V.value.

Returns:

Differentiable spike output with the same shape as input V. Values are approximately 0 (no spike) or 1 (spike) with smooth gradients for backpropagation.

Return type:

ArrayLike

Notes

The voltage is scaled to unitless values (mV) before applying the surrogate function. For Hodgkin-Huxley neurons, the actual spike threshold is V_T + 30 mV (default: -33 mV), but the surrogate function operates on the raw scaled voltage for gradient computation.

This method is primarily used for surrogate gradient learning. The discrete spike detection logic in the update method is independent and uses the three-condition test (refractory, threshold, local maximum).

Examples

>>> import saiunit as u
>>> import jax.numpy as jnp
>>> from brainpy_state import hh_cond_exp_traub
>>>
>>> neurons = hh_cond_exp_traub(10)
>>> neurons.init_state()
>>>
>>> # Get spike output for current state
>>> spikes = neurons.get_spike()
>>> print(spikes.shape)
(10,)
>>>
>>> # Get spike output for custom voltage
>>> V_custom = jnp.array([-60., -50., -40.]) * u.mV
>>> neurons_3 = hh_cond_exp_traub(3)
>>> neurons_3.init_state()
>>> spikes_custom = neurons_3.get_spike(V_custom)

See also

update

Main update method with discrete spike detection logic.

init_state(**kwargs)

Initialize all state variables for the neuron population.

Initializes membrane potential, gating variables, synaptic conductances, stimulation current buffer, refractory counter, and last spike time. If initial values are not explicitly provided, they are computed as follows:

• V: defaults to E_L

• m, h, n: computed from equilibrium at initial V using Traub-Miles rate equations (without V_T offset, matching NEST initialization)

• g_ex, g_in: initialized to zero

• I_stim: initialized to zero

• refractory_step_count: initialized to zero (not refractory)

• last_spike_time: initialized to -1e7 ms (far in the past)

Parameters:

**kwargs – Unused compatibility parameters accepted by the base-state API.

Notes

The equilibrium gating variable computation uses the raw voltage V (not V - V_T) to match NEST’s initialization procedure. During dynamics, the rate equations use the shifted voltage V - V_T, but initialization uses the unshifted value for consistency with NEST’s State_::State_ constructor.

This initialization ensures the neuron starts in a stable resting state when V_m_init = E_L (default). For custom initial voltages, gating variables are automatically adjusted to the corresponding equilibrium.

Examples

>>> import brainstate as bst
>>> import saiunit as u
>>> from brainpy_state import hh_cond_exp_traub
>>>
>>> # Initialize with default rest state
>>> neurons = hh_cond_exp_traub(100)
>>> neurons.init_state()
>>> print(neurons.V.value[0])  # Should be E_L = -60 mV
-60.0 mV
>>>
>>> # Initialize with custom voltage
>>> neurons = hh_cond_exp_traub(100, V_m_init=-65*u.mV)
>>> neurons.init_state()
>>> print(neurons.V.value[0])
-65.0 mV

Raises:

• ValueError – If an initializer cannot be broadcast to requested shape.

• TypeError – If initializer outputs have incompatible units/dtypes for the corresponding state variables.

See also

_hh_cond_exp_traub_equilibrium

Computes equilibrium gating values.

update(x=Quantity(0., 'pA'))

Update neuron state for one simulation step.

Integrates the 6-dimensional ODE system for one time step using adaptive RKF45 solver, processes incoming synaptic inputs, detects spikes based on threshold crossing and local maximum, and updates refractory state.

The update follows the NEST hh_cond_exp_traub update order:

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

2. Integrate the full 6-dimensional ODE system over one time step using an adaptive RKF45 solver.

3. Add arriving synaptic conductance jumps to g_ex / g_in.

4. Check spike condition: V_m >= V_T + 30 and V_old > V_m (threshold + local maximum).

5. Update refractory counter and record spike time.

6. Store buffered stimulation current for the next step.

Parameters:

x (ArrayLike, default 0 pA) – External stimulation current input (in addition to I_e), shape () or (*in_size,). This current is added to the constant I_e parameter and any registered current inputs via add_current_input().

Returns:

Spike output with shape (*in_size,). Values are computed using the surrogate spike function for differentiability. Spikes occur only when the discrete spike condition is satisfied (not refractory, threshold crossed, and local maximum detected).

Return type:

ArrayLike

Notes

Integration Details:

Each neuron’s state is integrated using an adaptive RKF45 integrator implemented in JAX with unit-aware arithmetic. This matches NEST’s GSL RKF45 solver. The ODE system is:

\[\begin{split}\frac{d}{dt}\begin{bmatrix} V_m \\ m \\ h \\ n \\ g_{ex} \\ g_{in} \end{bmatrix} = \begin{bmatrix} (-I_{Na} - I_K - I_L - I_{syn,ex} - I_{syn,in} + I_{stim} + I_e) / C_m \\ \alpha_m - (\alpha_m + \beta_m) m \\ \alpha_h - (\alpha_h + \beta_h) h \\ \alpha_n - (\alpha_n + \beta_n) n \\ -g_{ex} / \tau_{syn,ex} \\ -g_{in} / \tau_{syn,in} \end{bmatrix}\end{split}\]

Spike Detection Logic:

A spike is detected when all three conditions are met:

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

2. V_m >= V_T + 30 (threshold crossing)

3. V_old > V_m (local maximum - voltage falling)

No voltage reset occurs; repolarization is handled by intrinsic currents.

Synaptic Input Processing:

Delta inputs (spike events) are collected and split by sign:

• Positive weights -> excitatory conductance (g_ex += w)

• Negative weights -> inhibitory conductance (g_in += |w|)

Conductance jumps are applied after ODE integration, matching NEST’s update sequence.

Computational Complexity

Integration is performed with an adaptive vectorized RKF45 loop, including in-loop spike detection and refractory handling. All arithmetic is unit-aware via saiunit.math.

Failure Modes

• If the integrator detects numerical instability (V < -1e3 mV or V > 1e3 mV), a runtime error is raised.

• Extreme parameter values (very large conductances, very small time constants) may cause numerical instability.

See also

init_state

Initialize state variables.

get_spike

Compute surrogate spike output.