hh_psc_alpha_gap

class brainpy.state.hh_psc_alpha_gap(in_size, E_L=Quantity(-70., 'mV'), C_m=Quantity(40., 'pF'), g_Na=Quantity(4500., 'nS'), g_Kv1=Quantity(9., 'nS'), g_Kv3=Quantity(9000., 'nS'), g_L=Quantity(10., 'nS'), E_Na=Quantity(74., 'mV'), E_K=Quantity(-90., '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=None, Act_m_init=None, Inact_h_init=None, Act_n_init=None, Inact_p_init=None, spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', gsl_error_tol=1e-06, name=None)

NEST-compatible Hodgkin-Huxley neuron with alpha PSCs and gap junctions.

Short Description

Conductance-based spiking neuron model implementing Hodgkin-Huxley dynamics with alpha-function postsynaptic currents and support for gap-junction coupling. Uses modified ion-channel kinetics from Mancilla et al. (2007) with two distinct potassium conductances (Kv1, Kv3) for modeling gap-junction-coupled inhibitory interneurons.

Model Overview

hh_psc_alpha_gap extends the classic Hodgkin-Huxley formalism with:

• Sodium (Na) conductance: Activation gate \(m\), inactivation gate \(h\)

• Two potassium conductances: Fast Kv3 with \(p\) gate, slow Kv1 with \(n\) gate

• Leak conductance: Passive membrane current

• Alpha-function PSCs: Second-order synaptic current dynamics for excitatory and inhibitory inputs

• Gap-junction support: External resistive coupling current \(I_{gap}\)

• Hybrid spike detection: Combines voltage threshold (0 mV) with local maximum detection

• Explicit refractoriness: Suppresses spike emission during refractory period; subthreshold dynamics continue evolving

This implementation replicates NEST’s hh_psc_alpha_gap model (models/hh_psc_alpha_gap.{h,cpp}), using adaptive Runge-Kutta integration (RK45/Dormand-Prince) to match NEST’s GSL RKF45 solver.

1. Membrane Potential Dynamics

The membrane voltage evolves according to:

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

where the ionic currents are:

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

The potassium current combines contributions from slow Kv1 channels (\(n^4\) gating) and fast Kv3 channels (\(p^2\) gating), which is the key difference from standard HH models.

2. Gating Variable Dynamics

All four gating variables \(x \in \{m, h, n, p\}\) follow first-order kinetics:

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

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

\[\begin{split}\alpha_m &= \frac{40\,(V - 75.5)}{1 - e^{-(V - 75.5)/13.5}}, \quad \beta_m = \frac{1.2262}{e^{V/42.248}} \\ \alpha_h &= \frac{0.0035}{e^{V/24.186}}, \quad \beta_h = \frac{0.017\,(51.25 + V)}{1 - e^{-(51.25 + V)/5.2}} \\ \alpha_n &= \frac{0.014\,(V + 44)}{1 - e^{-(V + 44)/2.3}}, \quad \beta_n = \frac{0.0043}{e^{(V + 44)/34}} \\ \alpha_p &= \frac{V - 95}{1 - e^{-(V - 95)/11.8}}, \quad \beta_p = \frac{0.025}{e^{V/22.222}}\end{split}\]

These kinetics differ from the classic Hodgkin-Huxley equations and are based on experimental measurements from neocortical interneurons.

3. Gap-Junction Current

Gap junctions provide resistive electrical coupling:

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

where \(g_{ij}\) is the gap-junction conductance between neuron i and neuron j, and \(V_j\) is the membrane potential of the coupled neuron. In this single-neuron model, \(I_{gap}\) must be computed externally and provided as input via the x parameter to update() or through add_current_input().

4. Alpha-Function Synaptic Currents

Each synapse type (excitatory/inhibitory) uses a second-order system to generate alpha-shaped postsynaptic currents:

\[\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}\]

An incoming spike with weight \(w\) (in pA) increments \(dI_{syn}\) by \(w \cdot e / \tau_{syn}\), ensuring the peak current reaches \(w\) pA. The factor \(e = \exp(1)\) normalizes the alpha function.

5. Spike Detection Mechanism

Spikes are detected using a combined threshold-and-local-maximum criterion:

1. Not in refractory period: r == 0

2. Threshold crossing: \(V_m \geq 0\) mV

3. Local maximum: \(V_{old} > V_m\) (voltage is decreasing)

All three conditions must be satisfied simultaneously. This prevents multiple spike detections during the rising and falling phases of the action potential. Unlike integrate-and-fire models, no voltage reset occurs—repolarization happens naturally through activation of potassium currents.

6. Refractory Period

During the refractory period (duration \(t_{ref}\)), spike emission is suppressed, but the neuron’s subthreshold dynamics continue to evolve according to the differential equations. This differs from models that clamp the membrane potential during refractoriness.

7. Numerical Integration

NEST uses GSL’s RKF45 (Runge-Kutta-Fehlberg 4th/5th order) adaptive integrator with absolute tolerance 1e-6 and relative tolerance 0. This implementation uses a vectorized adaptive RKF45 integrator via AdaptiveRungeKuttaStep with matching tolerances. The 9-dimensional ODE system (V, m, h, n, p, dI_ex, I_ex, dI_in, I_in) is integrated simultaneously for all neurons over each time step.

Computational Complexity

• Per neuron, per time step: One adaptive ODE integration (~10-50 function evaluations depending on step size control)

• Scaling: Linear in population size (vectorized across neurons)

• Memory: O(population_size) for state storage

• JIT-compiled: Uses JAX-based adaptive RKF45 integrator for high performance and GPU compatibility

Parameters:

• in_size (int or tuple of int) – Population shape. Can be an integer (1D population) or tuple of integers (multidimensional population). Defines the number of neurons in the group.

• E_L (ArrayLike, optional) – Leak reversal potential (resting potential). Scalar or array with shape broadcastable to in_size. Unit: mV. Default: -70.0 mV.

• C_m (ArrayLike, optional) – Membrane capacitance. Must be strictly positive. Scalar or array with shape broadcastable to in_size. Unit: pF. Default: 40.0 pF.

• g_Na (ArrayLike, optional) – Sodium peak conductance. Must be non-negative. Scalar or array with shape broadcastable to in_size. Unit: nS. Default: 4500.0 nS.

• g_Kv1 (ArrayLike, optional) – Potassium Kv1 (slow) peak conductance. Must be non-negative. Scalar or array with shape broadcastable to in_size. Unit: nS. Default: 9.0 nS.

• g_Kv3 (ArrayLike, optional) – Potassium Kv3 (fast) peak conductance. Must be non-negative. Scalar or array with shape broadcastable to in_size. Unit: nS. Default: 9000.0 nS.

• g_L (ArrayLike, optional) – Leak conductance. Must be non-negative. Scalar or array with shape broadcastable to in_size. Unit: nS. Default: 10.0 nS.

• E_Na (ArrayLike, optional) – Sodium reversal potential. Scalar or array with shape broadcastable to in_size. Unit: mV. Default: 74.0 mV.

• E_K (ArrayLike, optional) – Potassium reversal potential. Scalar or array with shape broadcastable to in_size. Unit: mV. Default: -90.0 mV.

• t_ref (ArrayLike, optional) – Duration of refractory period. Must be non-negative. During this period, spike emission is suppressed but dynamics continue evolving. Scalar or array with shape broadcastable to in_size. Unit: ms. Default: 2.0 ms.

• tau_syn_ex (ArrayLike, optional) – Excitatory synaptic time constant (alpha-function rise time). Must be strictly positive. Scalar or array with shape broadcastable to in_size. Unit: ms. Default: 0.2 ms.

• tau_syn_in (ArrayLike, optional) – Inhibitory synaptic time constant (alpha-function rise time). Must be strictly positive. Scalar or array with shape broadcastable to in_size. Unit: ms. Default: 2.0 ms.

• I_e (ArrayLike, optional) – Constant external input current. Positive values are depolarizing. Scalar or array with shape broadcastable to in_size. Unit: pA. Default: 0.0 pA.

• V_m_init (ArrayLike or None, optional) – Initial membrane potential. If None, uses NEST’s default value of -69.604012 mV. Scalar or array with shape broadcastable to in_size. Unit: mV. Default: None.

• Act_m_init (ArrayLike or None, optional) – Initial sodium activation gating variable. Must be in [0, 1]. If None, computed from equilibrium at V_m_init. Scalar or array with shape broadcastable to in_size. Unitless. Default: None.

• Inact_h_init (ArrayLike or None, optional) – Initial sodium inactivation gating variable. Must be in [0, 1]. If None, computed from equilibrium at V_m_init. Scalar or array with shape broadcastable to in_size. Unitless. Default: None.

• Act_n_init (ArrayLike or None, optional) – Initial Kv1 activation gating variable. Must be in [0, 1]. If None, computed from equilibrium at V_m_init. Scalar or array with shape broadcastable to in_size. Unitless. Default: None.

• Inact_p_init (ArrayLike or None, optional) – Initial Kv3 activation gating variable. Must be in [0, 1]. If None, computed from equilibrium at V_m_init. Scalar or array with shape broadcastable to in_size. Unitless. Default: None.

• spk_fun (Callable, optional) – Surrogate gradient function for differentiable spike generation. Should be a callable from braintools.surrogate with signature (ArrayLike) -> ArrayLike. Used for gradient-based learning. Default: braintools.surrogate.ReluGrad().

• spk_reset ({'hard', 'soft'}, optional) – Spike reset mode. For HH models, this affects surrogate gradient computation only (no actual voltage reset occurs). ‘hard’: stop gradient propagation; ‘soft’: allow gradient flow. Default: ‘hard’.

• gsl_error_tol (float, optional) – Absolute tolerance for the embedded RKF45 error estimate. Must be strictly positive. Default: 1e-6 (matching NEST).

• name (str or None, optional) – Name of the neuron population for identification. If None, an automatic name is generated. Default: None.

Parameter Mapping

Parameter	Default	Math Symbol	Description
in_size	(required)	—	Population shape
E_L	-70.0 mV	\(E_L\)	Leak reversal potential (resting potential)
C_m	40.0 pF	\(C_m\)	Membrane capacitance
g_Na	4500.0 nS	\(g_{Na}\)	Sodium peak conductance
g_Kv1	9.0 nS	\(g_{Kv1}\)	Potassium Kv1 (slow) peak conductance
g_Kv3	9000.0 nS	\(g_{Kv3}\)	Potassium Kv3 (fast) peak conductance
g_L	10.0 nS	\(g_L\)	Leak conductance
E_Na	74.0 mV	\(E_{Na}\)	Sodium reversal potential
E_K	-90.0 mV	\(E_K\)	Potassium reversal potential
t_ref	2.0 ms	\(t_{ref}\)	Duration of refractory period
tau_syn_ex	0.2 ms	\(\tau_{syn,ex}\)	Excitatory synaptic time constant
tau_syn_in	2.0 ms	\(\tau_{syn,in}\)	Inhibitory synaptic time constant
I_e	0.0 pA	\(I_e\)	Constant external input current
V_m_init	-69.60401 mV	—	Initial membrane potential (NEST default)
Act_m_init	None	—	Initial Na activation (None -> equilibrium)
Inact_h_init	None	—	Initial Na inactivation (None -> equilibrium)
Act_n_init	None	—	Initial Kv1 activation (None -> equilibrium)
Inact_p_init	None	—	Initial Kv3 activation (None -> equilibrium)
spk_fun	ReluGrad()	—	Surrogate spike function
spk_reset	‘hard’	—	Reset mode for gradient computation
gsl_error_tol	1e-6	—	Absolute tolerance for RKF45 error estimate

V

Membrane potential \(V_m\). Shape: in_size. Unit: mV.

Type:

brainstate.HiddenState

m

Sodium activation gating variable. Shape: in_size. Range: [0, 1]. Unitless.

Type:

brainstate.HiddenState

h

Sodium inactivation gating variable. Shape: in_size. Range: [0, 1]. Unitless.

Type:

brainstate.HiddenState

n

Potassium Kv1 activation gating variable. Shape: in_size. Range: [0, 1]. Unitless.

Type:

brainstate.HiddenState

p

Potassium Kv3 activation gating variable. Shape: in_size. Range: [0, 1]. Unitless.

Type:

brainstate.HiddenState

I_syn_ex

Excitatory postsynaptic current. Shape: in_size. Unit: pA.

Type:

brainstate.ShortTermState

I_syn_in

Inhibitory postsynaptic current. Shape: in_size. Unit: pA.

Type:

brainstate.ShortTermState

dI_syn_ex

Excitatory alpha-kernel derivative state. Shape: in_size. Unit: pA/ms.

Type:

brainstate.ShortTermState

dI_syn_in

Inhibitory alpha-kernel derivative state. Shape: in_size. Unit: pA/ms.

Type:

brainstate.ShortTermState

I_stim

Stimulation current buffer for next time step. Shape: in_size. Unit: pA.

Type:

brainstate.ShortTermState

refractory_step_count

Refractory countdown in discrete time steps. Counts down from ceil(t_ref / dt) to 0. Shape: in_size. Unit: steps (integer).

Type:

brainstate.ShortTermState

integration_step

Persistent RKF45 substep size estimate (ms).

Type:

brainstate.ShortTermState

last_spike_time

Time of most recent spike. Shape: in_size. Unit: ms.

Type:

brainstate.ShortTermState

Raises:

• ValueError – If C_m <= 0 (capacitance must be strictly positive).

• ValueError – If t_ref < 0 (refractory time cannot be negative).

• ValueError – If tau_syn_ex <= 0 or tau_syn_in <= 0 (time constants must be strictly positive).

• ValueError – If any conductance (g_Na, g_Kv1, g_Kv3, g_L) is negative.

Notes

Differences from ``hh_psc_alpha``:

• Adds gap-junction current \(I_{gap}\) to membrane equation

• Uses modified ion-channel kinetics (Mancilla et al. 2007) with two potassium channel types (Kv1, Kv3)

• Different default conductance values optimized for interneuron models

Spike Weights and Input Interpretation:

• Positive spike weights -> excitatory input (added to dI_syn_ex)

• Negative spike weights -> inhibitory input (added to dI_syn_in)

• Weight magnitude in pA determines peak current amplitude

• Gap-junction current should be provided via x parameter to update() or registered via add_current_input()

Integration Accuracy:

• Adaptive step-size control ensures high accuracy but variable computational cost per step

• Default tolerance matches NEST for comparable numerical behavior

Gradient-Based Learning:

• Surrogate gradients enable backpropagation through spike generation

• The spk_fun parameter controls the shape of the surrogate gradient

• No actual voltage reset occurs, so gradients flow through natural action potential dynamics

References

[1]

Hodgkin AL, Huxley AF (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology 117:500-544. DOI: https://doi.org/10.1113/jphysiol.1952.sp004764

[2]

Mancilla JG, Lewis TG, Pinto DJ, Rinzel J, Connors BW (2007). Synchronization of electrically coupled pairs of inhibitory interneurons in neocortex. Journal of Neuroscience, 27:2058-2073. DOI: https://doi.org/10.1523/JNEUROSCI.2715-06.2007

[3]

Gerstner W, Kistler W (2002). Spiking neuron models: Single neurons, populations, plasticity. Cambridge University Press.

[4]

Hahne J, Helias M, Kunkel S, Igarashi J, Bolten M, Frommer A, Diesmann M (2015). A unified framework for spiking and gap-junction interactions in distributed neuronal network simulations. Frontiers in Neuroinformatics, 9:22. DOI: https://doi.org/10.3389/fninf.2015.00022

See also

hh_psc_alpha

Hodgkin-Huxley neuron without gap-junction support.

hh_cond_exp_traub

Alternative HH implementation with exponential PSCs.

iaf_cond_exp

Simpler integrate-and-fire model with conductance-based synapses.

Examples

Create a single gap-junction-coupled HH neuron:

>>> import brainpy.state as bs
>>> import saiunit as u
>>> neuron = bs.hh_psc_alpha_gap(in_size=1, E_L=-70*u.mV, C_m=40*u.pF)
>>> neuron.init_all_states()

Simulate with constant input current:

>>> import brainstate as bst
>>> with bst.environ.context(dt=0.1*u.ms):
...     neuron.init_all_states()
...     spikes = []
...     for i in range(1000):
...         spk = neuron.update(x=500*u.pA)  # 500 pA input
...         spikes.append(spk.item())

Create a population with heterogeneous capacitance:

>>> import jax.numpy as jnp
>>> C_m_values = jnp.linspace(30, 50, 10) * u.pF
>>> neurons = bs.hh_psc_alpha_gap(in_size=10, C_m=C_m_values)
>>> neurons.init_all_states()

Add gap-junction coupling between two neurons:

>>> neuron1 = bs.hh_psc_alpha_gap(in_size=1)
>>> neuron2 = bs.hh_psc_alpha_gap(in_size=1)
>>> neuron1.init_all_states()
>>> neuron2.init_all_states()
>>> g_gap = 0.5 * u.nS  # gap-junction conductance
>>> # In update loop:
>>> I_gap_1 = g_gap * (neuron2.V.value - neuron1.V.value)
>>> I_gap_2 = g_gap * (neuron1.V.value - neuron2.V.value)
>>> spk1 = neuron1.update(x=I_gap_1)
>>> spk2 = neuron2.update(x=I_gap_2)

get_spike(V=None)

Compute differentiable spike output from membrane potential.

Applies the surrogate gradient function (spk_fun) to generate a continuous, differentiable spike signal from the membrane potential. This enables gradient-based learning through the spiking dynamics.

For the HH model with combined threshold-and-local-maximum spike detection, this method is typically called with a specially crafted voltage signal (positive for spiking, negative otherwise) rather than the raw membrane potential.

Parameters:

V (ArrayLike or None, optional) – Membrane potential or voltage-derived signal. If None, uses the current value of self.V.value. Shape: in_size. Unit: mV. Default: None.

Returns:

spike_signal – Differentiable spike output with shape in_size. The surrogate function maps the scaled voltage to a continuous output (typically in range [0, 1] or [-1, 1] depending on spk_fun). Gradients flow through this function during backpropagation.

Return type:

ArrayLike

Notes

• The voltage is normalized by dividing by 1 mV before applying spk_fun to ensure dimensionless input

• For binary spike detection, threshold the returned value at 0

• The choice of spk_fun affects gradient magnitudes and learning dynamics (e.g., ReluGrad, SigmoidGrad, SuperSpike)

See also

update

Main simulation step that calls this method.

braintools.surrogate

Module containing surrogate gradient functions.

Examples

Direct spike computation:

>>> neuron = bs.hh_psc_alpha_gap(in_size=1)
>>> neuron.init_all_states()
>>> neuron.V.value = 10 * u.mV  # depolarized
>>> spk = neuron.get_spike()
>>> print(f"Spike signal: {spk}")

Using custom voltage signal:

>>> V_custom = jnp.array([1e-12, -1.0]) * u.mV  # spike/no-spike
>>> spk = neuron.get_spike(V=V_custom)

init_state(**kwargs)

Initialize all state variables for the neuron population.

Sets up membrane potential, gating variables, synaptic currents, and internal state tracking. Gating variables are initialized to their equilibrium values at the initial membrane potential unless explicitly specified, ensuring the neuron starts in a consistent resting state without transient artifacts.

Initialization Strategy:

1. Membrane potential: Set to V_m_init (default: NEST’s -69.604012 mV equilibrium value)

2. Gating variables: If Act_m_init, Inact_h_init, Act_n_init, or Inact_p_init are None, compute equilibrium values at V_m_init using _hh_psc_alpha_gap_equilibrium()

3. Synaptic currents: Initialize I_syn_ex, I_syn_in and their derivatives to zero

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

5. Spike timing: Set last_spike_time to large negative value (-1e7 ms)

Equilibrium Initialization Rationale:

Starting gating variables at their equilibrium values for the given initial voltage prevents spurious transient currents during the first few time steps. Without this, arbitrary initial values would cause artificial spikes or oscillations as the system relaxes to equilibrium.

Parameters:

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

Attributes Set

This method initializes the following instance attributes:

Vbrainstate.HiddenState

Membrane potential. Initial value: V_m_init. Unit: mV.

mbrainstate.HiddenState

Sodium activation gating variable. Initial value: Act_m_init if provided, otherwise equilibrium at V_m_init. Range: [0, 1].

hbrainstate.HiddenState

Sodium inactivation gating variable. Initial value: Inact_h_init if provided, otherwise equilibrium at V_m_init. Range: [0, 1].

nbrainstate.HiddenState

Kv1 activation gating variable. Initial value: Act_n_init if provided, otherwise equilibrium at V_m_init. Range: [0, 1].

pbrainstate.HiddenState

Kv3 activation gating variable. Initial value: Inact_p_init if provided, otherwise equilibrium at V_m_init. Range: [0, 1].

I_syn_exbrainstate.ShortTermState

Excitatory synaptic current. Initial value: 0. Unit: pA.

I_syn_inbrainstate.ShortTermState

Inhibitory synaptic current. Initial value: 0. Unit: pA.

dI_syn_exbrainstate.ShortTermState

Time derivative of excitatory current (alpha kernel state). Initial value: 0. Unit: pA/ms.

dI_syn_inbrainstate.ShortTermState

Time derivative of inhibitory current (alpha kernel state). Initial value: 0. Unit: pA/ms.

I_stimbrainstate.ShortTermState

Stimulation current buffer. Initial value: 0. Unit: pA.

refractory_step_countbrainstate.ShortTermState

Refractory countdown in time steps. Initial value: 0 (not refractory). Unit: steps (int32).

integration_stepbrainstate.ShortTermState

Persistent RKF45 substep size estimate (ms).

last_spike_timebrainstate.ShortTermState

Time of last spike. Initial value: -1e7 ms (far in the past). Unit: ms.

Notes

• Must be called before update() to ensure state variables exist

• Can be called multiple times to reinitialize (e.g., between trials)

• For heterogeneous populations with per-neuron initial conditions, pass arrays to V_m_init, Act_m_init, etc. during construction

• The NEST default initial voltage (-69.604012 mV) places the neuron near its resting state with minimal initial transients

State Variable Types:

• HiddenState: For slow variables (V, gating variables) that persist across time steps and require gradient tracking

• ShortTermState: For fast variables (currents, counters) that are recomputed each step or have short-term dynamics

See also

_hh_psc_alpha_gap_equilibrium

Computes equilibrium gating values.

update

Main simulation step that uses these state variables.

Examples

Basic initialization:

>>> neuron = bs.hh_psc_alpha_gap(in_size=10)
>>> neuron.init_state()
>>> print(neuron.V.value.shape)
(10,)

Custom initial conditions:

>>> import jax.numpy as jnp
>>> V_init = jnp.linspace(-75, -65, 10) * u.mV
>>> neuron = bs.hh_psc_alpha_gap(in_size=10, V_m_init=V_init)
>>> neuron.init_state()
>>> print(neuron.V.value)  # voltage varies across population

Initialize with custom gating variables:

>>> neuron = bs.hh_psc_alpha_gap(
...     in_size=1,
...     Act_m_init=0.1,  # specific sodium activation
...     Inact_h_init=0.9  # specific sodium inactivation
... )
>>> neuron.init_state()
>>> print(f"m={neuron.m.value}, h={neuron.h.value}")

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

Advance neuron state by one simulation time step.

Executes the full update cycle following NEST’s hh_psc_alpha_gap implementation order. This includes integrating the 9-dimensional ODE system (membrane potential, gating variables, and synaptic currents), processing incoming spikes, detecting output spikes, and managing the refractory state.

Update Sequence (Matching NEST Order):

1. Record pre-integration voltage: Save V_old for spike detection

2. Integrate ODEs: Solve the 9D system over [t, t+dt] using adaptive RKF45

3. Process arriving spikes: Add weighted spike inputs to dI_syn_ex / dI_syn_in derivative states

4. Detect spikes: Check three-part condition (not refractory, threshold crossed, local maximum)

5. Update refractory counter: Reset to t_ref/dt steps if spiking, otherwise decrement

6. Store stimulation buffer: Save I_stim for next time step

7. Return spike output: Compute surrogate spike signal

The 9-Dimensional ODE System:

The ODE integrator solves for state vector \(\mathbf{y} = [V_m, m, h, n, p, dI_{ex}, I_{ex}, dI_{in}, I_{in}]\) using the dynamics described in the class docstring. All neurons in the population are integrated simultaneously (vectorized).

Spike Input Processing:

• Spike inputs arrive via sum_delta_inputs() (collects all registered delta-function inputs)

• Positive weights -> excitatory: added to dI_syn_ex

• Negative weights -> inhibitory: added to dI_syn_in

• Normalization factor \(e/\tau_{syn}\) ensures peak current equals weight magnitude

Current Input Processing:

• Continuous current inputs via sum_current_inputs() (collects parameter x and all registered current inputs)

• Gap-junction current typically provided through x parameter

• Also includes constant bias current I_e

Spike Detection Logic:

spike = (r == 0) & (V_m >= 0.0) & (V_old > V_m)

This ensures only one spike per action potential by requiring: (1) not refractory, (2) above threshold, (3) voltage decreasing (local maximum has passed).

Numerical Considerations

• All neurons are integrated simultaneously via vectorized adaptive RKF45 (JIT-compiled, GPU-compatible)

• Adaptive step-size control may use 10-50 function evaluations per time step depending on dynamics and error tolerances

Integration Tolerances:

The ODE solver uses gsl_error_tol as absolute tolerance to control step-size adaptation. Smaller values increase accuracy but require more function evaluations. Default value (1e-6) matches NEST’s GSL settings.

Parameters:

x (ArrayLike, optional) – External input current. Can be scalar or array broadcastable to population shape. Typically includes gap-junction current computed as \(\sum_j g_{ij}(V_j - V_i)\) for coupled networks. Also accepts stimulation currents from external devices. Unit: pA. Default: 0.0 pA.

Returns:

spike_output – Differentiable spike signal with shape in_size. Computed by applying surrogate gradient function spk_fun to a voltage-derived signal: positive when spiking (V_out = 1e-12), negative otherwise (V_out = -1.0). For binary spike detection, threshold at 0. For gradient-based learning, use the returned analog values with spk_fun’s surrogate gradient.

Return type:

ArrayLike

Notes

State Update Side Effects:

This method modifies the following instance attributes:

• V.value: Updated membrane potential

• m.value, h.value, n.value, p.value: Updated gating variables

• I_syn_ex.value, I_syn_in.value: Updated synaptic currents

• dI_syn_ex.value, dI_syn_in.value: Updated synaptic derivatives

• I_stim.value: Buffered stimulation for next step

• refractory_step_count.value: Updated refractory countdown

• last_spike_time.value: Spike time when spiking occurs

Gap-Junction Usage Example:

For a network with gap-junction coupling matrix G and voltage vector V:

>>> G = [[0, 0.5], [0.5, 0]] * u.nS  # coupling conductances
>>> V1, V2 = neuron1.V.value, neuron2.V.value
>>> I_gap1 = G[0,1] * (V2 - V1)
>>> I_gap2 = G[1,0] * (V1 - V2)
>>> spk1 = neuron1.update(x=I_gap1)
>>> spk2 = neuron2.update(x=I_gap2)

Alternative Input Mechanism:

Instead of passing gap-junction current via x, you can register it as a named current input:

>>> neuron.add_current_input('gap', lambda: I_gap)
>>> spk = neuron.update()  # gap current applied automatically

Warnings

• Do not call update() before init_state() or init_all_states() — state variables must be initialized first

• Ensure time step dt is sufficiently small (typically <= 0.1 ms) for accurate spike detection and alpha-function dynamics

See also

init_state

Initialize neuron state variables.

get_spike

Compute spike output from membrane potential.

sum_delta_inputs

Collect all registered delta-function inputs.

sum_current_inputs

Collect all registered current inputs.