ht_neuron

class brainpy.state.ht_neuron(in_size, E_Na=30.0, E_K=-90.0, g_NaL=0.2, g_KL=1.0, tau_m=16.0, theta_eq=-51.0, tau_theta=2.0, tau_spike=1.75, t_ref=2.0, g_peak_AMPA=0.1, tau_rise_AMPA=0.5, tau_decay_AMPA=2.4, E_rev_AMPA=0.0, g_peak_NMDA=0.075, tau_rise_NMDA=4.0, tau_decay_NMDA=40.0, E_rev_NMDA=0.0, V_act_NMDA=-25.57, S_act_NMDA=0.081, tau_Mg_slow_NMDA=22.7, tau_Mg_fast_NMDA=0.68, instant_unblock_NMDA=False, g_peak_GABA_A=0.33, tau_rise_GABA_A=1.0, tau_decay_GABA_A=7.0, E_rev_GABA_A=-70.0, g_peak_GABA_B=0.0132, tau_rise_GABA_B=60.0, tau_decay_GABA_B=200.0, E_rev_GABA_B=-90.0, g_peak_NaP=1.0, E_rev_NaP=30.0, N_NaP=3.0, g_peak_KNa=1.0, E_rev_KNa=-90.0, tau_D_KNa=1250.0, g_peak_T=1.0, E_rev_T=0.0, N_T=2.0, g_peak_h=1.0, E_rev_h=-40.0, voltage_clamp=False, gsl_error_tol=0.001, spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', name=None)

NEST-compatible Hill-Tononi thalamocortical neuron model with intrinsic currents.

Implements the conductance-based integrate-and-fire neuron model from Hill & Tononi (2005) designed to simulate sleep-wake dynamics in thalamocortical networks. Features adaptive threshold, repolarizing post-spike potassium current, four receptor types (AMPA, NMDA, GABA_A, GABA_B), and four intrinsic currents (I_NaP, I_KNa, I_T, I_h) that mediate burst firing, adaptation, and oscillatory behavior.

This implementation replicates NEST’s ht_neuron (models/ht_neuron.{h,cpp}) using JAX-compatible adaptive ODE integration with AdaptiveRungeKuttaStep.

Parameters:

• in_size (int or tuple of int) – Population shape (e.g., 100 or (10, 10)). Determines the number of neurons in this layer.

• E_Na (float, default 30.0) – Sodium reversal potential in mV. Sets the depolarized reset level after spike.

• E_K (float, default -90.0) – Potassium reversal potential in mV. Sets the hyperpolarized target for repolarization during refractory period.

• g_NaL (float, default 0.2) – Sodium leak conductance (unitless). Contributes depolarizing leak current.

• g_KL (float, default 1.0) – Potassium leak conductance (unitless). Contributes hyperpolarizing leak current.

• tau_m (float, default 16.0) – Membrane time constant in ms. Governs the rate of subthreshold voltage changes.

• theta_eq (float, default -51.0) – Equilibrium spike threshold in mV. The threshold relaxes to this value with time constant tau_theta.

• tau_theta (float, default 2.0) – Threshold relaxation time constant in ms. Controls adaptation timescale.

• tau_spike (float, default 1.75) – Repolarization time constant for post-spike potassium current in ms. Governs the speed of voltage recovery during refractory period.

• t_ref (float, default 2.0) – Absolute refractory period in ms. Duration of post-spike potassium current.

• g_peak_AMPA (float, default 0.1) – Peak AMPA conductance (unitless). Scaled by spike inputs to produce excitatory synaptic current.

• tau_rise_AMPA (float, default 0.5) – AMPA conductance rise time constant in ms. Must be < tau_decay_AMPA.

• tau_decay_AMPA (float, default 2.4) – AMPA conductance decay time constant in ms.

• E_rev_AMPA (float, default 0.0) – AMPA reversal potential in mV.

• g_peak_NMDA (float, default 0.075) – Peak NMDA conductance (unitless). Subject to voltage-dependent Mg²⁺ block.

• tau_rise_NMDA (float, default 4.0) – NMDA conductance rise time constant in ms. Must be < tau_decay_NMDA.

• tau_decay_NMDA (float, default 40.0) – NMDA conductance decay time constant in ms.

• E_rev_NMDA (float, default 0.0) – NMDA reversal potential in mV.

• V_act_NMDA (float, default -25.57) – Voltage at 50% NMDA Mg²⁺ unblock in mV.

• S_act_NMDA (float, default 0.081) – Slope of NMDA Mg²⁺ unblocking sigmoid in 1/mV.

• tau_Mg_slow_NMDA (float, default 22.7) – Slow Mg²⁺ unblocking time constant in ms. Must be > tau_Mg_fast_NMDA.

• tau_Mg_fast_NMDA (float, default 0.68) – Fast Mg²⁺ unblocking time constant in ms.

• instant_unblock_NMDA (bool, default False) – If True, use instantaneous Mg²⁺ unblocking (m^NMDA = m_∞). If False, use two-stage kinetics with fast and slow unblocking components.

• g_peak_GABA_A (float, default 0.33) – Peak GABA_A conductance (unitless). Fast inhibitory synaptic current.

• tau_rise_GABA_A (float, default 1.0) – GABA_A rise time constant in ms. Must be < tau_decay_GABA_A.

• tau_decay_GABA_A (float, default 7.0) – GABA_A decay time constant in ms.

• E_rev_GABA_A (float, default -70.0) – GABA_A reversal potential in mV.

• g_peak_GABA_B (float, default 0.0132) – Peak GABA_B conductance (unitless). Slow inhibitory synaptic current.

• tau_rise_GABA_B (float, default 60.0) – GABA_B rise time constant in ms. Must be < tau_decay_GABA_B.

• tau_decay_GABA_B (float, default 200.0) – GABA_B decay time constant in ms.

• E_rev_GABA_B (float, default -90.0) – GABA_B reversal potential in mV.

• g_peak_NaP (float, default 1.0) – Peak persistent sodium current conductance (unitless). Mediates subthreshold depolarization and bistability.

• E_rev_NaP (float, default 30.0) – I_NaP reversal potential in mV.

• N_NaP (float, default 3.0) – I_NaP activation exponent (power to which m_∞ is raised).

• g_peak_KNa (float, default 1.0) – Peak I_KNa conductance (unitless). Provides slow spike-frequency adaptation.

• E_rev_KNa (float, default -90.0) – I_KNa reversal potential in mV.

• tau_D_KNa (float, default 1250.0) – I_KNa D-variable relaxation time constant in ms. Large value produces very slow adaptation (~seconds).

• g_peak_T (float, default 1.0) – Peak low-threshold Ca²⁺ current conductance (unitless). Mediates rebound bursts and oscillations.

• E_rev_T (float, default 0.0) – I_T reversal potential in mV.

• N_T (float, default 2.0) – I_T activation exponent (power to which m_IT is raised).

• g_peak_h (float, default 1.0) – Peak hyperpolarization-activated current conductance (unitless). Contributes to rebound excitation and resonance.

• E_rev_h (float, default -40.0) – I_h reversal potential in mV.

• voltage_clamp (bool, default False) – If True, clamp membrane potential at its initial value throughout simulation. Used for testing intrinsic current dynamics in isolation.

• gsl_error_tol (float, default 1e-3) – Absolute error tolerance for the adaptive RKF45 integrator.

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

• spk_reset (str, default 'hard') – Spike reset mode: ‘hard’ (stop gradient) or ‘soft’ (V -= V_th).

• name (str or None, default None) – Name of this neuron population.

Parameter Mapping

The table below maps brainpy.state parameter names to NEST equivalents:

brainpy.state Parameter	NEST Parameter	Default	Units
in_size	(N/A)	—	—
E_Na	E_Na	30.0	mV
E_K	E_K	-90.0	mV
g_NaL	g_NaL	0.2	(unitless)
g_KL	g_KL	1.0	(unitless)
tau_m	tau_m	16.0	ms
theta_eq	theta_eq	-51.0	mV
tau_theta	tau_theta	2.0	ms
tau_spike	tau_spike	1.75	ms
t_ref	t_ref	2.0	ms
g_peak_AMPA	g_peak_AMPA	0.1	(unitless)
tau_rise_AMPA	tau_rise_AMPA	0.5	ms
tau_decay_AMPA	tau_decay_AMPA	2.4	ms
E_rev_AMPA	E_rev_AMPA	0.0	mV
g_peak_NMDA	g_peak_NMDA	0.075	(unitless)
tau_rise_NMDA	tau_rise_NMDA	4.0	ms
tau_decay_NMDA	tau_decay_NMDA	40.0	ms
E_rev_NMDA	E_rev_NMDA	0.0	mV
V_act_NMDA	V_act_NMDA	-25.57	mV
S_act_NMDA	S_act_NMDA	0.081	1/mV
tau_Mg_slow_NMDA	tau_Mg_slow_NMDA	22.7	ms
tau_Mg_fast_NMDA	tau_Mg_fast_NMDA	0.68	ms
instant_unblock_NMDA	instant_unblock	False	—
g_peak_GABA_A	g_peak_GABA_A	0.33	(unitless)
tau_rise_GABA_A	tau_rise_GABA_A	1.0	ms
tau_decay_GABA_A	tau_decay_GABA_A	7.0	ms
E_rev_GABA_A	E_rev_GABA_A	-70.0	mV
g_peak_GABA_B	g_peak_GABA_B	0.0132	(unitless)
tau_rise_GABA_B	tau_rise_GABA_B	60.0	ms
tau_decay_GABA_B	tau_decay_GABA_B	200.0	ms
E_rev_GABA_B	E_rev_GABA_B	-90.0	mV
g_peak_NaP	g_peak_NaP	1.0	(unitless)
E_rev_NaP	E_rev_NaP	30.0	mV
N_NaP	NaP_N	3.0	—
g_peak_KNa	g_peak_KNa	1.0	(unitless)
E_rev_KNa	E_rev_KNa	-90.0	mV
tau_D_KNa	tau_D_KNa	1250.0	ms
g_peak_T	g_peak_T	1.0	(unitless)
E_rev_T	E_rev_T	0.0	mV
N_T	T_N	2.0	—
g_peak_h	g_peak_h	1.0	(unitless)
E_rev_h	E_rev_h	-40.0	mV
voltage_clamp	voltage_clamp	False	—
gsl_error_tol	(GSL tolerance)	1e-3	—

Notes

1. Model Architecture

The ht_neuron is an integrate-and-fire model with:

• Adaptive threshold: Threshold increases transiently after spike, then relaxes to theta_eq, providing spike-frequency adaptation on ~ms timescale.

• Soft reset: No hard voltage reset. Instead, V and theta are set to E_Na, and a repolarizing K⁺ current drives voltage back toward E_K during t_ref.

• Four synaptic receptor types: AMPA (fast excitation), NMDA (slow excitation with voltage-dependent Mg²⁺ block), GABA_A (fast inhibition), GABA_B (slow inhibition). Each uses beta-function (biexponential) conductance time course.

• Four intrinsic currents:

  • I_NaP (persistent Na⁺): Subthreshold depolarizing current; enables bistability and up-states.

  • I_KNa (depolarization-activated K⁺): Very slow adaptation (~1 s timescale).

  • I_T (low-threshold Ca²⁺): Mediates rebound bursts; deinactivates during hyperpolarization and activates rapidly on depolarization.

  • I_h (hyperpolarization-activated cation current): Sag current; contributes to rebound and resonance.

2. Membrane Dynamics

The membrane potential obeys:

\[\frac{dV}{dt} = \frac{I_{leak} + I_{syn} + I_{intrinsic} + I_{stim}}{\tau_m} + I_{spike}\]

where:

\[\begin{split}I_{leak} &= -g_{NaL}(V - E_{Na}) - g_{KL}(V - E_K) \\ I_{syn} &= -g_{AMPA}(V - E_{AMPA}) - g_{NMDA} m^{NMDA}(V - E_{NMDA}) \\ &\quad - g_{GABA_A}(V - E_{GABA_A}) - g_{GABA_B}(V - E_{GABA_B}) \\ I_{intrinsic} &= I_{NaP} + I_{KNa} + I_T + I_h \\ I_{spike} &= \begin{cases} -(V - E_K) / \tau_{spike} & \text{if refractory} \\ 0 & \text{otherwise} \end{cases}\end{split}\]

3. Dynamic Threshold

\[\frac{d\theta}{dt} = -\frac{\theta - \theta_{eq}}{\tau_\theta}\]

On spike, theta is reset to E_Na (along with V), then decays back to theta_eq. This provides fast spike-frequency adaptation.

4. Beta-Function Synapses

Each synapse type uses a two-variable beta function (difference of exponentials):

\[\begin{split}\frac{dg'}{dt} &= -\frac{g'}{\tau_{rise}} \\ \frac{dg}{dt} &= g' - \frac{g}{\tau_{decay}}\end{split}\]

On arrival of a spike, the DG variable (g’) is incremented by:

\[\Delta g' = g_{peak} \cdot \text{norm}(\tau_{rise}, \tau_{decay}) \cdot w\]

where norm is the beta normalization factor and w is the synaptic weight.

5. NMDA Voltage Dependence

NMDA channels are blocked by Mg²⁺ at hyperpolarized potentials and unblock with depolarization. Two modes are available:

• Instantaneous unblocking (instant_unblock_NMDA=True):

  \[m^{NMDA} = \frac{1}{1 + \exp(-S_{act}(V - V_{act}))}\]

• Two-stage kinetics (instant_unblock_NMDA=False):

  \[\begin{split}\frac{dm_{fast}}{dt} &= \frac{m_\infty - m_{fast}}{\tau_{Mg,fast}} \\ \frac{dm_{slow}}{dt} &= \frac{m_\infty - m_{slow}}{\tau_{Mg,slow}} \\ m^{NMDA} &= A_1(V) m_{fast} + A_2(V) m_{slow}\end{split}\]

  where A₁(V) = 0.51 - 0.0028V and A₂ = 1 - A₁. This captures the experimentally observed slow Mg²⁺ unblocking kinetics (Vargas-Caballero & Robinson, 2003).

6. Intrinsic Currents

• I_NaP (persistent sodium):

  \[\begin{split}m_\infty &= \frac{1}{1 + \exp(-(V + 55.7)/7.7)} \\ I_{NaP} &= -g_{NaP} \cdot (m_\infty)^{N_{NaP}} \cdot (V - E_{NaP})\end{split}\]

  No inactivation; provides tonic depolarizing drive.

• I_KNa (depolarization-activated potassium):

  \[\begin{split}D_{influx} &= \frac{0.025}{1 + \exp(-(V + 10)/5)} \\ \frac{dD}{dt} &= \frac{\tau_{D,KNa} \cdot D_{influx} + 0.001 - D}{\tau_{D,KNa}} \\ m_\infty &= \frac{1}{1 + (0.25/D)^{3.5}} \\ I_{KNa} &= -g_{KNa} \cdot m_\infty \cdot (V - E_{KNa})\end{split}\]

  D accumulates slowly during depolarization; provides adaptation on ~second timescale.

• I_T (low-threshold Ca²⁺):

  \[\begin{split}m_\infty &= \frac{1}{1 + \exp(-(V + 59)/6.2)} \\ h_\infty &= \frac{1}{1 + \exp((V + 83)/4)} \\ \tau_m &= 0.22 / (\exp(-(V+132)/16.7) + \exp((V+16.8)/18.2)) + 0.13 \\ \tau_h &= 8.2 + \frac{56.6 + 0.27 \exp((V+115.2)/5)}{1 + \exp((V+86)/3.2)} \\ I_T &= -g_T \cdot m^{N_T} \cdot h \cdot (V - E_T)\end{split}\]

  Activation is fast; inactivation is slower. Channel deinactivates during hyperpolarization, enabling rebound bursts.

• I_h (hyperpolarization-activated cation current):

  \[\begin{split}m_\infty &= \frac{1}{1 + \exp((V + 75)/5.5)} \\ \tau_m &= \frac{1}{\exp(-14.59 - 0.086V) + \exp(-1.87 + 0.0701V)} \\ I_h &= -g_h \cdot m \cdot (V - E_h)\end{split}\]

  Activates slowly at hyperpolarized potentials; provides depolarizing sag and contributes to rebound.

7. Spike Detection and Reset

A spike occurs when ref_steps == 0 and V >= theta. On spike:

• V → E_Na (≈ +30 mV)

• theta → E_Na

• ref_steps → ceil(t_ref / dt) + 1

During the refractory period, I_spike drives V back toward E_K.

8. Numerical Integration

The model uses AdaptiveRungeKuttaStep with RKF45 (Runge-Kutta-Fehlberg 4(5) adaptive integration). This matches NEST’s GSL RKF45 integrator in terms of order and adaptive step-size control.

9. Conductance Units

All conductances are unitless in this model. The membrane equation is written as dV/dt = I/tau_m, meaning currents have units of mV/ms (i.e., they are already divided by capacitance). Peak conductances g_peak_* scale the synaptic currents.

10. Sleep-Wake Transitions

The ht_neuron was designed to model thalamocortical neurons that exhibit two distinct firing modes:

• Tonic firing (awake/depolarized): Regular spiking driven by excitatory input and I_NaP.

• Burst firing (sleep/hyperpolarized): Rebound bursts mediated by I_T deinactivation and I_h rebound.

By varying the balance of intrinsic current conductances (g_peak_T, g_peak_h, g_peak_NaP) and background synaptic input, the model can transition between these modes, reproducing the sleep-wake dynamics observed in thalamocortical circuits.

Examples

Example 1: Single neuron with injected current

>>> import brainpy as bp
>>> import brainpy.state as bps
>>> import saiunit as u
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>>
>>> # Create a single ht_neuron
>>> neuron = bps.ht_neuron(1, g_peak_T=0.0, g_peak_h=0.0)
>>>
>>> # Initialize state
>>> with bp.environ.context(dt=0.1 * u.ms):
...     neuron.init_all_states()
...
...     # Simulate 500 ms with step current
...     currents = np.concatenate([
...         np.zeros(1000),
...         np.ones(3000) * 2.0,  # 2 mV/ms injected current
...         np.zeros(1000)
...     ])
...     voltages = []
...     for I in currents:
...         neuron.update(I)
...         voltages.append(neuron.V.value[0])
>>>
>>> # Plot membrane potential
>>> times = np.arange(len(voltages)) * 0.1
>>> plt.figure(figsize=(10, 4))
>>> plt.plot(times, voltages)
>>> plt.xlabel('Time (ms)')
>>> plt.ylabel('Membrane potential (mV)')
>>> plt.title('ht_neuron response to step current')
>>> plt.show()

Example 2: Rebound burst with I_T

>>> # Enable I_T for burst firing
>>> neuron = bps.ht_neuron(1, g_peak_T=1.0, g_peak_h=0.5)
>>>
>>> with bp.environ.context(dt=0.1 * u.ms):
...     neuron.init_all_states()
...
...     # Hyperpolarize, then release
...     currents = np.concatenate([
...         np.zeros(500),
...         -np.ones(1000) * 3.0,  # hyperpolarizing current
...         np.zeros(1500)
...     ])
...     voltages = []
...     for I in currents:
...         neuron.update(I)
...         voltages.append(neuron.V.value[0])
>>>
>>> # Observe rebound burst after hyperpolarization ends
>>> plt.figure(figsize=(10, 4))
>>> plt.plot(np.arange(len(voltages)) * 0.1, voltages)
>>> plt.xlabel('Time (ms)')
>>> plt.ylabel('V (mV)')
>>> plt.title('Rebound burst mediated by I_T')
>>> plt.show()

Example 3: Multi-receptor synaptic input

>>> # Network with AMPA and NMDA receptors
>>> pre = bps.LIF(100, V_rest=-70*u.mV, V_th=-50*u.mV, V_reset=-70*u.mV,
...               tau=20*u.ms, R=1*u.ohm, V_initializer=bp.init.Normal(-70, 5))
>>> post = bps.ht_neuron(10, g_peak_AMPA=0.1, g_peak_NMDA=0.05,
...                      instant_unblock_NMDA=False)
>>>
>>> # Create projections for AMPA and NMDA
>>> ampa_proj = bps.AlignPostProj(
...     pre=pre, post=post,
...     comm=bp.event.FixedProb(0.1, weight=1.0),
...     syn=bps.Expon.desc(tau=2.4 * u.ms),
...     label='AMPA'
... )
>>> nmda_proj = bps.AlignPostProj(
...     pre=pre, post=post,
...     comm=bp.event.FixedProb(0.1, weight=1.0),
...     syn=bps.Expon.desc(tau=40 * u.ms),
...     label='NMDA'
... )
>>>
>>> # Simulate network dynamics
>>> # (implementation depends on BrainPy network API)

See also

hh_psc_alpha

Hodgkin-Huxley neuron with alpha-shaped PSCs

iaf_cond_exp

Simple IAF with exponential conductance synapses

aeif_cond_alpha

Adaptive exponential IAF with alpha conductances

References

[1]

Hill S, Tononi G (2005). Modeling sleep and wakefulness in the thalamocortical system. Journal of Neurophysiology, 93:1671-1698. DOI: https://doi.org/10.1152/jn.00915.2004

[2]

Vargas-Caballero M, Robinson HPC (2003). A slow fraction of Mg²⁺ unblock of NMDA receptors limits their contribution to spike generation in cortical pyramidal neurons. Journal of Neurophysiology, 89:2778-2783. DOI: https://doi.org/10.1152/jn.01038.2002

get_spike(V=None)

Generate differentiable spike output using surrogate gradient function.

Converts the discrete spike condition (V >= theta) into a continuous, differentiable output suitable for gradient-based optimization. The voltage is scaled relative to the dynamic threshold before passing through the surrogate function.

Parameters:

V (ArrayLike or None, default None) – Membrane potential in mV. If None, uses self.V.value. Can be a scalar, 1D, or multi-dimensional array matching the neuron population shape. For explicit spike injection (e.g., after reset), pass a manually set value (e.g., large positive for spike, large negative for no spike).

Returns:

Surrogate spike output with the same shape as V. Values near 1.0 indicate spike, values near 0.0 indicate no spike. Gradients flow through the surrogate function (e.g., ReluGrad, sigmoid, etc.) rather than being zero.

Return type:

ArrayLike

Notes

1. Voltage Scaling

The input voltage is scaled by the threshold magnitude to normalize the surrogate function input:

\[v_{scaled} = \frac{V - \theta}{\max(|\theta_{eq}|, 1)}\]

This ensures that v_scaled ≈ 0 when V ≈ theta, and v_scaled > 0 when spiking. The denominator prevents numerical issues if theta_eq is very small.

2. Surrogate Function

The scaled voltage is passed through the surrogate gradient function specified during initialization (default: braintools.surrogate.ReluGrad()):

\[s = \text{spk\_fun}(v_{scaled})\]

During forward pass, this typically produces a Heaviside-like step (0 or 1). During backward pass, the gradient is replaced by a smooth approximation (e.g., d/dv max(0, v) = 1 if v > 0 else 0 for ReluGrad).

3. Spike Detection vs. Spike Output

This method generates the output for surrogate gradient learning. The actual spike detection (threshold crossing, reset, refractory logic) happens in the update() method using discrete logic. The two are synchronized: when update() detects a spike, it calls get_spike() with a manually set V value to ensure the output is 1.0.

4. Gradient Flow

Unlike a true Heaviside function (which has zero gradient everywhere except at the discontinuity), the surrogate function provides non-zero gradients in a neighborhood around the threshold, enabling backpropagation through spiking networks.

Examples

>>> import brainpy.state as bps
>>> import saiunit as u
>>> import brainstate
>>>
>>> neuron = bps.ht_neuron(1)
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     neuron.init_all_states()
...
...     # Check spike output at rest (V ≈ -70 mV, theta ≈ -51 mV)
...     # V < theta, so no spike
...     spike = neuron.get_spike()
...     print(spike)  # ≈ 0.0
...
...     # Manually set V above threshold
...     neuron.V.value = -45.0  # > theta
...     spike = neuron.get_spike()
...     print(spike)  # ≈ 1.0 (depends on surrogate function)

init_state(**kwargs)

Initialize all state variables to physiologically consistent equilibrium values.

Sets the membrane potential to the leak reversal potential (weighted average of E_Na and E_K based on leak conductances), threshold to theta_eq, all synaptic variables to zero, and all intrinsic gating variables to their voltage-dependent equilibrium values at the initial membrane potential.

Parameters:

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

Notes

1. Initial Membrane Potential

Computed from leak conductance balance:

\[V_{init} = \frac{g_{NaL} \cdot E_{Na} + g_{KL} \cdot E_K}{g_{NaL} + g_{KL}}\]

With default parameters (g_NaL=0.2, g_KL=1.0, E_Na=30 mV, E_K=-90 mV):

\[V_{init} = \frac{0.2 \cdot 30 + 1.0 \cdot (-90)}{0.2 + 1.0} = -70\ \text{mV}\]

2. Threshold Initialization

theta is set to theta_eq (default -51 mV).

3. Synaptic Variables

All beta-function state variables are initialized to zero:

• DG_AMPA, G_AMPA = 0

• DG_NMDA, G_NMDA = 0

• DG_GABA_A, G_GABA_A = 0

• DG_GABA_B, G_GABA_B = 0

4. Intrinsic Gating Variables

All gating variables are set to their steady-state values at V_init:

• m_fast_NMDA = m_slow_NMDA = m_∞^NMDA(V_init)

• m_Ih = m_∞^Ih(V_init)

• D_IKNa = D_∞(V_init)

• m_IT = m_∞^IT(V_init)

• h_IT = h_∞^IT(V_init)

At V_init ≈ -70 mV (resting potential):

• m_Ih ≈ 0.4 (partially activated, since I_h activates at hyperpolarization)

• m_IT ≈ 0.05 (mostly deactivated)

• h_IT ≈ 0.9 (mostly deinactivated, ready to support burst)

• D_IKNa ≈ 0.001 (minimal adaptation at rest)

• m_NMDA ≈ 0.01 (strongly blocked by Mg²⁺ at rest)

5. Refractory Counter

ref_steps = 0 (neuron is not refractory at initialization).

6. Stimulation Current

I_stim = 0 (no external current at t=0).

7. Spike Time

last_spike_time = -1e7 ms (far in the past, ensures no artificial refractory period at simulation start).

8. Voltage Clamp Value

If voltage_clamp=True, _V_clamp is set to V_init and will be enforced during all subsequent updates.

This initialization matches NEST’s ht_neuron::State_::set() and calibrate() methods, ensuring consistent starting conditions for simulation comparisons.

update(x=0.0)

Advance neuron state by one simulation time step with adaptive ODE integration.

Performs a complete update cycle for the ht_neuron model, including: (1) adaptive RKF45 integration of the 16-dimensional ODE system, (2) spike detection and reset, (3) refractory period management, (4) synaptic input processing, and (5) external current buffering. The update sequence matches NEST’s ht_neuron::update() implementation for numerical consistency.

Parameters:

x (float or ArrayLike, default 0.0) –

External stimulation current in mV/ms (since conductances are unitless in this model, currents are expressed as I/C_m). Can be:

• Scalar: Applied uniformly to all neurons in the population.

• Array: Shape must broadcast to the neuron population shape (e.g., for spatially varying input or per-neuron stimulation protocols).

This current is added to the membrane equation as I_stim and affects dV/dt during the next integration step.

Returns:

Differentiable spike output with shape matching the neuron population. Values near 1.0 indicate spike, near 0.0 indicate no spike. Compatible with surrogate gradient-based learning.

Return type:

ArrayLike

Notes

1. Update Sequence

The update follows this precise ordering (matching NEST):

Step 1: ODE Integration

Integrate the 16-dimensional state vector from t to t+dt using adaptive RKF45 via AdaptiveRungeKuttaStep. The state vector contains:

\[\mathbf{y} = [V, \theta, DG_{AMPA}, G_{AMPA}, DG_{NMDA}, G_{NMDA}, DG_{GABA_A}, G_{GABA_A}, DG_{GABA_B}, G_{GABA_B}, m_{fast}^{NMDA}, m_{slow}^{NMDA}, m_{Ih}, D_{IKNa}, m_{IT}, h_{IT}]\]

The ODE right-hand side (defined in _vector_field) computes:

• Membrane potential derivative from leak, synaptic, intrinsic, and stimulation currents, plus post-spike repolarization if refractory

• Threshold relaxation: dθ/dt = -(θ - θ_eq)/tau_θ

• Beta-function synaptic conductance dynamics (4 receptor types)

• NMDA Mg²⁺ unblocking kinetics (fast and slow components)

• Intrinsic gating variable dynamics (I_h, I_T, I_KNa)

Step 2: Post-Integration Constraints

After integration, enforce:

• Voltage clamp: If voltage_clamp=True, reset V to _V_clamp.

• Instantaneous NMDA blocking: Clamp m_fast_NMDA and m_slow_NMDA to not exceed m_∞^NMDA(V), ensuring the Mg²⁺ block cannot be “overshot” during adaptive time steps.

Step 3: Spike Detection and Reset

If ref_steps == 0 and V >= theta, a spike is generated:

• V → E_Na (≈ +30 mV)

• θ → E_Na

• ref_steps → ceil(t_ref / dt) + 1

• spike_flag = True

Step 4: Refractory Counter Decrement

If ref_steps > 0, decrement by 1. This happens after spike detection, so a neuron that just spiked will spend t_ref ms refractory.

Step 5: Synaptic Spike Input Delivery

Add arriving spikes to the DG (derivative of conductance) variables. Inputs are retrieved from delta_inputs with labels ‘AMPA’, ‘NMDA’, ‘GABA_A’, ‘GABA_B’:

\[DG_{receptor} \mathrel{+}= g_{peak,receptor} \cdot \text{norm} \cdot w \cdot N_{spikes}\]

Unlabeled delta inputs default to AMPA.

Step 6: Stimulation Current Buffering

Store the input current x in I_stim for use in the next update cycle. This matches NEST’s one-step delay for external currents.

2. Refractory Dynamics

During the refractory period (ref_steps > 0), the neuron cannot spike, and the post-spike potassium current is active:

\[I_{spike} = -\frac{V - E_K}{\tau_{spike}}\]

This drives V toward E_K (hyperpolarization) with time constant tau_spike.

3. Synaptic Input Routing

The ht_neuron expects delta inputs to be labeled by receptor type. Projections should add inputs via:

post.add_delta_input(weight * pre_spike, label='AMPA')

If no label is provided, inputs accumulate in the generic delta_inputs and are routed to AMPA by default.

4. Numerical Considerations

• Adaptive integration: The RKF45 solver uses variable step sizes to maintain accuracy. Typical internal steps are ~0.01-0.1 ms depending on voltage dynamics.

• Vectorized integration: All neurons in the population are integrated simultaneously using JAX vectorized operations via AdaptiveRungeKuttaStep.

• Intrinsic current caching: Intrinsic currents (I_NaP, I_KNa, I_T, I_h) are computed after integration and stored in separate state variables for recording.

5. Gradient Compatibility

The integration uses JAX-based AdaptiveRungeKuttaStep, enabling automatic differentiation through the integration. Combined with surrogate gradient spike output, the model supports end-to-end backpropagation.

Warnings

• Unlabeled inputs default to AMPA: If you send synaptic inputs without specifying a receptor label, they will be routed to AMPA receptors by default. This may produce unexpected results if you intended NMDA, GABA_A, or GABA_B.