hh_cond_beta_gap_traub#

class brainpy.state.hh_cond_beta_gap_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(-50., 'mV'), E_ex=Quantity(0., 'mV'), E_in=Quantity(-80., 'mV'), t_ref=Quantity(2., 'ms'), tau_rise_ex=Quantity(0.5, 'ms'), tau_decay_ex=Quantity(5., 'ms'), tau_rise_in=Quantity(0.5, 'ms'), tau_decay_in=Quantity(10., 'ms'), I_e=Quantity(0., 'pA'), gsl_error_tol=1e-06, V_m_init=None, Act_m_init=None, Inact_h_init=None, Act_n_init=None, spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', name=None)#

NEST-compatible Hodgkin-Huxley neuron with beta-function synapses and gap junctions.

Implements a conductance-based Hodgkin-Huxley model with Traub-Miles gating kinetics, beta-function (double-exponential) synaptic conductances, and support for gap-junction coupling. Based on the NEST hh_cond_beta_gap_traub model.

1. Model Overview

This model extends the classical Hodgkin-Huxley formalism to include:

  • Traub-Miles gating kinetics: Simplified three-variable (\(m\), \(h\), \(n\)) sodium and potassium channel dynamics from Traub and Miles (1991) [1].

  • Beta-function synapses: Double-exponential conductance profiles with separate rise and decay time constants for excitatory and inhibitory inputs.

  • Gap-junction support: Resistive coupling current that can be supplied externally to model electrical synapses between neurons.

  • Refractory spike detection: Physiological spike detection based on threshold crossing and local maximum detection with refractory period enforcement.

This is a point neuron model (single compartment) suitable for large-scale network simulations where detailed morphology is not required but synaptic dynamics and gap-junction coupling are important.

2. Membrane Dynamics

The membrane potential evolves according to:

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

where the ionic and synaptic 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}\]

Physical interpretation:

  • \(I_{Na}\) – Fast sodium current responsible for spike upstroke.

  • \(I_K\) – Delayed rectifier potassium current for repolarization.

  • \(I_L\) – Leak current maintaining resting potential.

  • \(I_{syn,ex}\), \(I_{syn,in}\) – Excitatory and inhibitory synaptic currents.

  • \(I_{gap}\) – Gap-junction current from electrically coupled neighbors.

3. Gating Variable Dynamics

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

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

with Traub-Miles rate functions using voltage-shifted dynamics \(V = V_m - V_T\):

\[\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 -50 mV) effectively shifts the spike threshold.

Computational note: Singularities in \(\alpha\) functions at specific voltages are handled via L’Hôpital’s rule in the ODE solver.

4. Beta-Function Synaptic Conductances

Synaptic conductances follow double-exponential (beta-function) dynamics:

\[\begin{split}\frac{d(\Delta g_{ex})}{dt} &= -\frac{\Delta g_{ex}}{\tau_{decay,ex}} \\ \frac{dg_{ex}}{dt} &= \Delta g_{ex} - \frac{g_{ex}}{\tau_{rise,ex}}\end{split}\]

and analogously for inhibitory conductance \(g_{in}\).

Spike input handling:

  • Excitatory spikes (positive weights) increment \(\Delta g_{ex}\).

  • Inhibitory spikes (negative weights) increment \(\Delta g_{in}\) (sign-flipped).

  • Each spike adds \(w \times \text{PSConInit}\) to \(\Delta g\), where \(\text{PSConInit}\) is the beta normalization factor ensuring peak conductance of 1 nS for unit weight.

Why beta functions? Unlike simple exponential or alpha functions, beta functions provide independent control over rise and decay time scales, critical for accurately modeling AMPA (fast), NMDA (slow), and GABA receptors.

5. Gap-Junction Current

Gap junctions model electrical synapses as resistive couplings:

\[I_{gap} = \sum_j g_{gap,ij}\,(V_j - V_i)\]

In this single-neuron implementation, \(I_{gap}\) must be computed externally (e.g., by a network simulation framework) and supplied via the x parameter to update() or via add_current_input.

6. Spike Detection

A spike is emitted when all three conditions are satisfied:

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

  2. \(V_m \geq V_T + 30\) mV (threshold crossing)

  3. \(V_{old} > V_m\) (local maximum detection)

No voltage reset occurs after spike emission (unlike integrate-and-fire models); repolarization is driven naturally by the potassium current.

Refractory period: During refractory steps, spike emission is suppressed but subthreshold dynamics continue to evolve. This prevents multiple spike detections during the falling phase of an action potential.

7. Numerical Integration

Uses an adaptive Runge-Kutta-Fehlberg (RKF45) integrator implemented in JAX. Default absolute tolerance (gsl_error_tol=1e-6) matches NEST’s GSL RKF45 integrator settings for numerical correspondence in benchmark comparisons.

The ODE system has 8 state variables per neuron: \([V_m, m, h, n, \Delta g_{ex}, g_{ex}, \Delta g_{in}, g_{in}]\).

Parameters:
  • in_size (Size) – Shape of the neuron population. Can be int (1D), tuple of ints (multidimensional), or None (scalar neuron). Determines state variable array dimensions.

  • E_L (ArrayLike, default -60 mV) – Leak reversal potential (resting potential in absence of input).

  • C_m (ArrayLike, default 200 pF) – Membrane capacitance. Must be strictly positive. Typical range: 50-500 pF.

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

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

  • g_L (ArrayLike, default 10 nS) – Leak conductance. Must be non-negative. Determines input resistance and time constant.

  • E_Na (ArrayLike, default 50 mV) – Sodium reversal potential. Typically +40 to +60 mV.

  • E_K (ArrayLike, default -90 mV) – Potassium reversal potential. Typically -80 to -100 mV.

  • V_T (ArrayLike, default -50 mV) – Voltage offset for gating dynamics. Shifts the effective spike threshold.

  • E_ex (ArrayLike, default 0 mV) – Excitatory synaptic reversal potential (typical for AMPA/NMDA receptors).

  • E_in (ArrayLike, default -80 mV) – Inhibitory synaptic reversal potential (typical for GABA receptors).

  • t_ref (ArrayLike, default 2 ms) – Refractory period duration. Must be non-negative. Increase if multiple spikes are detected per action potential.

  • tau_rise_ex (ArrayLike, default 0.5 ms) – Excitatory synaptic rise time constant. Must be strictly positive.

  • tau_decay_ex (ArrayLike, default 5.0 ms) – Excitatory synaptic decay time constant. Must be strictly positive. Should be larger than tau_rise_ex for physiological beta-function shape.

  • tau_rise_in (ArrayLike, default 0.5 ms) – Inhibitory synaptic rise time constant. Must be strictly positive.

  • tau_decay_in (ArrayLike, default 10.0 ms) – Inhibitory synaptic decay time constant. Must be strictly positive.

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

  • gsl_error_tol (ArrayLike, default 1e-6) – Unitless local RKF45 error tolerance, broadcastable and strictly positive.

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

  • Act_m_init (ArrayLike, optional) – Initial sodium activation variable. If None, computed from equilibrium at V_m_init.

  • Inact_h_init (ArrayLike, optional) – Initial sodium inactivation variable. If None, computed from equilibrium at V_m_init.

  • Act_n_init (ArrayLike, optional) – Initial potassium activation variable. If None, computed from equilibrium at V_m_init.

  • spk_fun (Callable, default braintools.surrogate.ReluGrad()) – Surrogate gradient function for differentiable spike generation during backpropagation. Only affects gradient computation; forward-pass spike detection is always threshold-based.

  • spk_reset (str, default 'hard') – Spike reset mode for gradient computation. 'hard' uses stop_gradient; 'soft' allows gradients through spike. Does not affect forward dynamics.

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

Parameter Mapping

This table maps brainpy.state parameter names to NEST hh_cond_beta_gap_traub parameter names and mathematical symbols:

brainpy.state

NEST

Math

Description

in_size

(population size)

Number/shape of neurons

E_L

E_L

\(E_L\)

Leak reversal potential (mV)

C_m

C_m

\(C_m\)

Membrane capacitance (pF)

g_Na

g_Na

\(g_{Na}\)

Sodium conductance (nS)

g_K

g_K

\(g_K\)

Potassium conductance (nS)

g_L

g_L

\(g_L\)

Leak conductance (nS)

E_Na

E_Na

\(E_{Na}\)

Sodium reversal (mV)

E_K

E_K

\(E_K\)

Potassium reversal (mV)

V_T

V_T

\(V_T\)

Voltage offset (mV)

E_ex

E_ex

\(E_{ex}\)

Excitatory reversal (mV)

E_in

E_in

\(E_{in}\)

Inhibitory reversal (mV)

t_ref

t_ref

\(t_{ref}\)

Refractory period (ms)

tau_rise_ex

tau_rise_ex

\(\tau_{rise,ex}\)

Excitatory rise time (ms)

tau_decay_ex

tau_decay_ex

\(\tau_{decay,ex}\)

Excitatory decay time (ms)

tau_rise_in

tau_rise_in

\(\tau_{rise,in}\)

Inhibitory rise time (ms)

tau_decay_in

tau_decay_in

\(\tau_{decay,in}\)

Inhibitory decay time (ms)

I_e

I_e

\(I_e\)

External current (pA)

gsl_error_tol

RKF45 local error tolerance

V_m_init

(initial V_m)

\(V_m(t=0)\)

Initial membrane potential (mV)

Act_m_init

(initial Act_m)

\(m(t=0)\)

Initial Na activation (0-1)

Inact_h_init

(initial Inact_h)

\(h(t=0)\)

Initial Na inactivation (0-1)

Act_n_init

(initial Act_n)

\(n(t=0)\)

Initial K activation (0-1)

V#

Membrane potential in mV. Shape: (*in_size,).

Type:

brainstate.HiddenState

m#

Sodium activation gating variable (unitless, 0-1 range).

Type:

brainstate.HiddenState

h#

Sodium inactivation gating variable (unitless, 0-1 range).

Type:

brainstate.HiddenState

n#

Potassium activation gating variable (unitless, 0-1 range).

Type:

brainstate.HiddenState

dg_ex#

Time derivative of excitatory conductance in nS/ms (beta-function intermediate state).

Type:

brainstate.ShortTermState

g_ex#

Excitatory synaptic conductance in nS.

Type:

brainstate.HiddenState

dg_in#

Time derivative of inhibitory conductance in nS/ms (beta-function intermediate state).

Type:

brainstate.ShortTermState

g_in#

Inhibitory synaptic conductance in nS.

Type:

brainstate.HiddenState

I_stim#

Buffered stimulation current in pA (applied in next time step).

Type:

brainstate.ShortTermState

refractory_step_count#

Integer countdown of remaining refractory steps. Zero means neuron can spike.

Type:

brainstate.ShortTermState

integration_step#

Persistent RKF45 substep size estimate (ms).

Type:

brainstate.ShortTermState

last_spike_time#

Time of most recent spike emission in ms (for recording/analysis).

Type:

brainstate.ShortTermState

Raises:

ValueError – If C_m <= 0, t_ref < 0, any time constant <= 0, or any conductance < 0.

Notes

Usage Considerations:

  1. Synaptic weight units: Spike weights are interpreted in conductance units (nS). A weight of 1.0 produces a peak conductance of 1 nS at the synapse’s rise time.

  2. Excitatory vs. inhibitory synapses: The sign of the synaptic weight determines the receptor type:

    • Positive weights drive g_ex (excitatory, reversal at E_ex).

    • Negative weights drive g_in (inhibitory, reversal at E_in).

    The sign is automatically handled by _sum_signed_delta_inputs().

  3. Gap-junction current: Must be computed externally and provided via the x parameter to update() or registered with add_current_input. In a network, compute as \(\sum_j g_{gap,ij}(V_j - V_i)\) where \(V_j\) are neighbor potentials and \(g_{gap,ij}\) are coupling conductances.

  4. No voltage reset: Unlike integrate-and-fire models, the membrane potential is not reset after spike emission. The potassium current naturally drives repolarization and hyperpolarization.

  5. Refractory period tuning: If the model emits multiple spikes per action potential, increase t_ref. Traub and Miles (1991) used 3 ms; NEST defaults to 2 ms.

  6. Numerical stability: The adaptive RKF45 integrator handles the stiff HH dynamics robustly. If you encounter instability, try reducing gsl_error_tol or increasing the simulation time step dt.

  7. Performance: All neurons are integrated in a single vectorized adaptive RKF45 loop via JAX, providing efficient GPU/TPU execution.

Examples

Basic single-neuron simulation with step current:

>>> import brainpy.state as bst
>>> import saiunit as u
>>> import brainstate
>>> import matplotlib.pyplot as plt
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     neuron = bst.hh_cond_beta_gap_traub(1)
...     neuron.init_all_states()
...     # Apply 500 pA step current for 100 ms
...     times, voltages = [], []
...     for t in range(1000):  # 100 ms simulation
...         if 200 <= t < 700:  # 20-70 ms
...             neuron.update(500 * u.pA)
...         else:
...             neuron.update(0 * u.pA)
...         times.append(brainstate.environ.get('t'))
...         voltages.append(neuron.V.value.item())
>>> plt.plot(times, voltages)
>>> plt.xlabel('Time (ms)')
>>> plt.ylabel('Membrane potential (mV)')
>>> plt.title('HH neuron with step current input')
>>> plt.show()

Network simulation with gap junctions:

>>> # Two coupled neurons with gap junction
>>> neuron_pop = bst.hh_cond_beta_gap_traub(2, I_e=200 * u.pA)
>>> neuron_pop.init_all_states()
>>> g_gap = 50.0  # nS gap conductance
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     for _ in range(1000):
...         V = neuron_pop.V.value
...         # Compute gap currents: I_gap[i] = g_gap * (V[j] - V[i])
...         I_gap = u.math.zeros_like(V)
...         I_gap = I_gap.at[0].set(g_gap * u.nS * (V[1] - V[0]))
...         I_gap = I_gap.at[1].set(g_gap * u.nS * (V[0] - V[1]))
...         neuron_pop.update(I_gap)

Beta-function synapse with different time constants:

>>> # Slow NMDA-like synapse (tau_rise=2ms, tau_decay=50ms)
>>> neuron = bst.hh_cond_beta_gap_traub(
...     1,
...     tau_rise_ex=2.0 * u.ms,
...     tau_decay_ex=50.0 * u.ms,
... )
>>> neuron.init_all_states()
>>> # Add excitatory spike input at t=10ms
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     for t in range(1000):
...         if t == 100:  # t=10ms
...             neuron.delta_inputs['spike'] = lambda: 5.0 * u.nS
...         neuron.update()

See also

hh_cond_exp_traub

Hodgkin-Huxley Traub model with single-exponential synapses.

hh_psc_alpha_gap

Hodgkin-Huxley model with gap junctions and alpha-function PSCs.

hh_psc_alpha

Classic HH model with current-based alpha-function synapses.

References

get_spike(V=None)[source]#

Compute differentiable spike output from membrane potential.

Applies the surrogate gradient function (spk_fun) to the membrane potential to generate a differentiable spike signal for gradient-based learning. This is used internally by update() to compute the return value.

Forward Pass vs. Backward Pass:

  • Forward pass: Returns a binary-like spike indicator (1.0 where spike occurred, 0.0 otherwise) based on the three-condition spike detection in update().

  • Backward pass: Gradients flow through the surrogate function (e.g., ReluGrad), which provides a smooth approximation of the Heaviside step function.

Why Surrogate Gradients?

The true spike detection logic (threshold + local maximum + refractory) is non-differentiable. Surrogate gradient methods replace the zero-everywhere gradient of the Heaviside function with a smooth proxy (e.g., ReLU, sigmoid, exponential) during backpropagation, enabling gradient-based optimization of spiking networks.

Parameters:

V (ArrayLike, optional) – Membrane potential in millivolts. If None, uses self.V.value (current state). Shape must match (*in_size,).

Returns:

spike – Differentiable spike output with same shape as V. Forward values are approximately binary (close to 0 or 1); backward gradients are provided by the surrogate function.

Return type:

ArrayLike

Notes

Voltage Scaling:

The membrane potential is divided by 1 mV to convert from physical units to a unitless scale before passing to spk_fun. This ensures the surrogate function operates on dimensionless voltage values (typically in the range -80 to +50 for biological neurons).

Surrogate Function Choice:

The default braintools.surrogate.ReluGrad() uses a rectified linear gradient:

\[\begin{split}\\text{forward}(V) &= H(V) \quad \\text{(Heaviside step function)} \\\\ \\frac{d}{dV}\\text{backward}(V) &= \\begin{cases} 1 & \\text{if } V > 0 \\\\ 0 & \\text{otherwise} \\end{cases}\end{split}\]

Other options include:

  • Sigmoid(): Smooth logistic gradient.

  • Gaussian(): Gaussian-shaped gradient.

  • PiecewiseQuadratic(): Quadratic spline gradient.

See braintools.surrogate for available functions.

Spike Reset Mode:

The spk_reset parameter ('hard' or 'soft') controls whether gradients flow through the spike in update():

  • 'hard': Uses jax.lax.stop_gradient to prevent gradients from propagating through the spike event. Gradient flow stops at the spike.

  • 'soft': Allows gradients to flow through the spike (no stop_gradient). This can help learning but may be less biologically plausible.

This method does not directly apply spk_reset; it is handled in update().

Examples

Direct spike computation from voltage:

>>> import brainpy.state as bst
>>> import saiunit as u
>>> neuron = bst.hh_cond_beta_gap_traub(1)
>>> neuron.init_all_states()
>>> # Manually set voltage above threshold
>>> V_test = (-50 + 30 + 1) * u.mV  # V_T + 30 + 1 = -19 mV
>>> spike = neuron.get_spike(V_test)
>>> print(f"Spike value: {spike.item():.3f}")

Using custom surrogate function:

>>> import braintools
>>> neuron = bst.hh_cond_beta_gap_traub(
...     1,
...     spk_fun=braintools.surrogate.Sigmoid(alpha=5.0),
... )
>>> neuron.init_all_states()
>>> spike = neuron.get_spike(neuron.V.value)

See also

update

Main update method that uses this function to compute spike output.

braintools.surrogate

Module containing surrogate gradient functions.

init_state(**kwargs)[source]#

Initialize all state variables to equilibrium or user-specified values.

Sets up hidden states (membrane potential, gating variables, synaptic conductances) and short-term states (refractory counter, spike time buffer). By default, initializes to physiologically realistic equilibrium values matching NEST’s initialization protocol.

Initialization Protocol:

  1. Membrane potential: Defaults to E_L (resting potential) if V_m_init is None.

  2. Gating variables: If Act_m_init, Inact_h_init, or Act_n_init are None, compute equilibrium values \(x_{\infty} = \alpha_x / (\alpha_x + \beta_x)\) at the initial membrane potential without V_T offset (matching NEST).

  3. Synaptic conductances: Initialize dg_ex, g_ex, dg_in, g_in to zero.

  4. Refractory state: Set refractory_step_count to 0 (not refractory).

  5. Spike time: Set last_spike_time to -1e7 ms (no recent spike).

Parameters:

**kwargs (dict, optional) – Unused compatibility parameters accepted by the base-state API.

Notes

Equilibrium Computation:

The equilibrium gating variables are computed using the Traub-Miles rate functions evaluated at the raw initial voltage (without V_T offset):

\[\begin{split}m_{\infty} &= \\frac{\\alpha_m(V_0)}{\\alpha_m(V_0) + \\beta_m(V_0)} \\\\ h_{\infty} &= \\frac{\\alpha_h(V_0)}{\\alpha_h(V_0) + \\beta_h(V_0)} \\\\ n_{\infty} &= \\frac{\\alpha_n(V_0)}{\\alpha_n(V_0) + \\beta_n(V_0)}\end{split}\]

where \(V_0 =\) V_m_init (or E_L if None). This matches NEST’s State_::State_(const Parameters_&) constructor, which uses y_[0] (= E_L) without applying the V_T shift used during dynamics integration.

Why No V_T Offset?

The V_T offset is applied during dynamics integration (in the ODE right-hand side) to shift the effective spike threshold. However, equilibrium initialization uses the absolute membrane potential to ensure consistency with the model’s resting state before any dynamics occur.

Custom Initialization:

To initialize with specific gating variable values (e.g., after depolarization):

>>> neuron = bst.hh_cond_beta_gap_traub(
...     10,
...     V_m_init=-50 * u.mV,      # Depolarized initial state
...     Act_m_init=0.3,            # Custom Na activation
...     Inact_h_init=0.4,          # Custom Na inactivation
...     Act_n_init=0.2,            # Custom K activation
... )
>>> neuron.init_all_states()

State Variable Summary:

After calling init_state, the following attributes are available:

  • V (HiddenState): Membrane potential (mV)

  • m (HiddenState): Sodium activation (0-1)

  • h (HiddenState): Sodium inactivation (0-1)

  • n (HiddenState): Potassium activation (0-1)

  • dg_ex (ShortTermState): Excitatory conductance derivative (nS/ms)

  • g_ex (HiddenState): Excitatory conductance (nS)

  • dg_in (ShortTermState): Inhibitory conductance derivative (nS/ms)

  • g_in (HiddenState): Inhibitory conductance (nS)

  • I_stim (ShortTermState): Stimulation current buffer (pA)

  • refractory_step_count (ShortTermState): Refractory countdown (int)

  • integration_step (ShortTermState): RKF45 substep size (ms)

  • last_spike_time (ShortTermState): Last spike time (ms)

Examples

Default equilibrium initialization:

>>> import brainpy.state as bst
>>> neuron = bst.hh_cond_beta_gap_traub(5)
>>> neuron.init_all_states()
>>> print(neuron.V.value)  # Should be E_L = -60 mV
>>> print(neuron.m.value)  # Equilibrium at -60 mV

Custom depolarized initial state:

>>> neuron = bst.hh_cond_beta_gap_traub(
...     1,
...     V_m_init=-45 * u.mV,  # Near threshold
... )
>>> neuron.init_all_states()
>>> print(neuron.V.value)  # -45 mV
>>> print(neuron.m.value)  # Equilibrium at -45 mV (higher than at -60 mV)
update(x=Quantity(0., 'pA'))[source]#

Update neuron state for one simulation time step.

Advances all state variables by one time step dt following the NEST hh_cond_beta_gap_traub update protocol. Integrates the 8D ODE system using adaptive RKF45, applies synaptic conductance jumps, detects spikes, and updates refractory state.

Update Protocol (Matching NEST Order):

  1. Record pre-integration voltage: Store \(V_{old} = V_m(t)\) for spike detection (local maximum criterion).

  2. ODE integration: Integrate the 8-variable system \([V_m, m, h, n, \Delta g_{ex}, g_{ex}, \Delta g_{in}, g_{in}]\) from \(t\) to \(t + dt\) using adaptive RKF45.

  3. Apply synaptic inputs: Add arriving spike-triggered conductance jumps:

    \[\begin{split}\Delta g_{ex} &\leftarrow \Delta g_{ex} + w_{ex} \times \text{PSConInit}_{ex} \\ \Delta g_{in} &\leftarrow \Delta g_{in} + w_{in} \times \text{PSConInit}_{in}\end{split}\]

    where \(\text{PSConInit}\) is the beta normalization factor ensuring peak conductance of 1 nS for unit weight.

  4. Spike detection: Emit spike if all conditions are met:

    • refractory_step_count == 0 (not refractory)

    • \(V_m(t+dt) \geq V_T + 30\) mV (threshold crossed)

    • \(V_{old} > V_m(t+dt)\) (local maximum detected)

  5. Refractory state update: If spike detected, set refractory_step_count to \(\lceil t_{ref} / dt \rceil\); otherwise decrement if nonzero.

  6. Buffer next stimulation current: Store I_stim for next step (one-step delay matching NEST buffer semantics).

Parameters:

x (ArrayLike, default 0 pA) –

External stimulation current for this time step. This is added to I_e and should include:

  • Gap-junction current: \(I_{gap} = \sum_j g_{gap,ij}(V_j - V_i)\)

  • Any additional bias or time-varying input current

Shape must broadcast with (*in_size,) or be scalar. Unit: picoamperes (pA).

Returns:

spike – Binary spike output with shape (*in_size,). Dtype: float64. Values of 1.0 indicate at least one spike event occurred during the integrated interval \((t, t+dt]\).

Return type:

ArrayLike

Notes

Numerical Integration Details:

  • All neurons are integrated simultaneously using a vectorized adaptive RKF45 loop implemented in JAX, providing efficient GPU/TPU execution.

  • The RKF45 (Runge-Kutta-Fehlberg) method uses adaptive step-size control with error tolerance gsl_error_tol.

  • Integration includes in-loop spike detection with local maximum criterion.

Spike Detection Logic:

The three-condition spike criterion prevents multiple detections per action potential:

  1. Refractory guard: Ensures minimum inter-spike interval.

  2. Threshold crossing: Voltage must exceed \(V_T + 30\) mV.

  3. Local maximum: \(V_{old} > V_m\) ensures detection only at peak, not during rising or falling phases.

This physiological detection method differs from IF models’ threshold-reset mechanism.

Gap-Junction Current Handling:

Gap junctions are not computed internally. You must:

  1. Compute neighbor voltage differences externally.

  2. Calculate \(I_{gap,i} = \sum_j g_{gap,ij}(V_j - V_i)\).

  3. Pass the result as the x parameter.

For networks, this typically requires gathering \(V_j\) from connected neurons before calling update().

See also

init_state

Initialize state variables before calling update().

get_spike

Compute differentiable spike output from voltage.