CUBANet

Article 1 derived the continuous-time leaky integrate-and-fire (LIF) neuron from Maxwell’s equations and the membrane circuit. This article picks up where that derivation leaves off and asks: how do you turn the continuous LIF into a discrete update rule that runs on a GPU? Two integrator choices (snnTorch-style and zero-order hold) give two rungs — standard-snn and cuba — both backed by a single Python class, CUBANet, in src/pinglab/models.py. Every empirical result in the notebooks on the CUBA / standard-snn family traces back to the math here.

The implementation lives in src/pinglab/models.py. This article keeps the derivations.

Architecture

CUBANet (current-based LIF) is a single class with one switch:

discretisation ∈ {snntorch, zoh} — how the continuous LIF is discretised. Spikes drive the membrane directly (no synaptic filter); the conductance-filtered story lives one rung up in COBANet.

The two canonical rungs:

rung	discretisation
standard-snn	snnTorch
cuba	ZOH

Tutorial mode (—kaiming-init) is a further switch that overrides Dale’s law and weight scaling to match snn.Leaky bit-for-bit at initialisation, used only by the parity probe.

The LIF primitive

The forward pass and backward pass share a primitive snn_lif_step that performs one timestep of decay + drive + spike + reset:

mem = β * mem + I
s   = σ_surr(mem − θ)
if reset == "subtract":
    mem = mem − θ · s
else:  # "zero"
    mem = mem · (1 − s)

In words: decay the old membrane, add the (already-scaled) drive $I$ , threshold the new membrane to emit a spike, and apply the reset. The caller is responsible for scaling $I$ — that’s where the two discretisations diverge.

snnTorch discretisation

Used by standard-snn. The update is the snnTorch convention — decay times the previous membrane plus the unscaled drive — with the reset applied after the spike is emitted:

U_{t+1} = \beta\, U_t + W\, s_t + b, \qquad S_{t+1} = \sigma_{\text{surr}}(U_{t+1} - \theta) \tag{5}

U_{t+1} \leftarrow \begin{cases} 0 & \text{reset=zero, } S_{t+1} = 1 \\ U_{t+1} - \theta\, S_{t+1} & \text{reset=subtract} \end{cases} \tag{6}

Two things to notice. First, the spike is computed on the post-update membrane $U_{t+1}$ , not the previous one — that’s the order of operations in snn_lif_step. Second, $W$ and $b$ enter the update without a $\Delta t$ scaling, so the same trained weights produce different per-step contributions at different $\Delta t$ . This is the dt-fragility shared with snn.Leaky.

Bias-balloon derivation

Take the expectation of (5) under a stationary mean input rate $r$ (Hz). Each step, the expected input is $\mathbb{E}[s_t] = r\Delta t$ , so $\mathbb{E}[W s_t] = W r \Delta t$ . The fixed-point of the membrane is the value $U_{ss}$ that satisfies $U_{ss} = \beta U_{ss} + W r \Delta t + b$ , giving

U_{ss} = \frac{W r \Delta t + b}{1-\beta} \tag{7}

For $\Delta t \ll \tau$ , use the first-order expansion $1 - \beta \approx \Delta t / \tau$ :

U_{ss} \approx W r \tau \;+\; \frac{b \tau}{\Delta t} \tag{8}

The spike contribution $W r \tau$ is $\Delta t$ -invariant: halve $\Delta t$ and the per-step kick halves, but you take twice as many steps. The bias contribution scales as $1/\Delta t$ : halve $\Delta t$ and the steady-state bias contribution doubles. This is why standard-snn with a non-zero bias drifts toward firing-rate saturation as $\Delta t$ shrinks at eval time.

ZOH discretisation

Used by cuba. Start from the continuous-time LIF in canonical form ( $E_L = 0$ , $R = 1$ ):

\tau\, \frac{dV}{dt} = -V + I(t), \qquad I(t) = \sum_k W\, \delta(t - t_k) + b \tag{9}

where the spike train $\sum_k \delta(t - t_k)$ represents incoming presynaptic spikes at times $\{t_k\}$ .

Step 1 — set up an integrating factor

Rewrite (9) as

\frac{dV}{dt} + \frac{1}{\tau}\, V = \frac{I(t)}{\tau}

This is a first-order linear ODE.

Aside — what is an integrating factor

For any first-order linear ODE of the form $\frac{dy}{dt} + p(t)\,y = q(t)$ , an integrating factor $\mu(t)$ is a function chosen specifically so that multiplying both sides by $\mu(t)$ turns the left-hand side into a single perfect derivative $\frac{d}{dt}[\mu(t)\,y(t)]$ . The benefit is concrete: once the LHS is a total derivative, both sides can be integrated directly — the LHS by inspection ( $\int \frac{d}{dt}[\mu y]\,dt = \mu y + \text{const}$ ), the RHS as a definite integral of a known function. The non-trivial ODE collapses to one antiderivative.

For the LHS to collapse this way, we need $\mu(t)$ to satisfy

\mu(t)\,\frac{dy}{dt} + \mu(t)\,p(t)\,y \;=\; \frac{d}{dt}\!\left[\mu(t)\,y(t)\right] \;=\; \mu(t)\,\frac{dy}{dt} + \mu'(t)\,y.

Matching the two RHS gives the defining condition $\mu'(t) = p(t)\,\mu(t)$ , a separable ODE with solution $\mu(t) = e^{\int p(t)\,dt}$ (up to an unimportant multiplicative constant). In our problem $p(t) = 1/\tau$ is constant, so $\int p\,dt = t/\tau$ and

\mu(t) \;=\; e^{t/\tau}.

This is the unique (up to scale) function whose logarithmic derivative is $1/\tau$ — i.e. the function whose growth rate matches the membrane’s leak rate.

Multiplying both sides of the original ODE by $\mu(t)$ :

e^{t/\tau}\, \frac{dV}{dt} + \frac{1}{\tau}\, e^{t/\tau}\, V = \frac{I(t)}{\tau}\, e^{t/\tau}

The left-hand side is, by the product rule, the derivative of $V(t)\, e^{t/\tau}$ :

\frac{d}{dt}\!\left[V(t)\, e^{t/\tau}\right] = \frac{d V}{dt}\, e^{t/\tau} + V(t)\, \frac{1}{\tau}\, e^{t/\tau}

So the ODE becomes

\frac{d}{dt}\!\left[V(t)\, e^{t/\tau}\right] = \frac{I(t)}{\tau}\, e^{t/\tau} \tag{10a}

Step 2 — integrate from $t$ to $t + \Delta t$

Apply $\int_t^{t+\Delta t} (\cdot)\, ds$ to both sides of (10a).

Aside — why we rename the variable from

t

s

Equation (10a) holds at every instant of time, with $t$ used both as the independent variable inside the derivative and as a generic label. To integrate over an interval we need to make a distinction: one symbol that ranges over the interval (the dummy / bound variable), and two fixed symbols that label the endpoints. If we kept everything called $t$ , expressions like $\int_t^{t+\Delta t} V(t)\, dt$ would be ambiguous — is the $t$ in $V(t)$ the moving one or the fixed lower endpoint?

The standard fix is to rename the dummy variable. The relation (10a) was

\frac{d}{dt}\!\left[V(t)\, e^{t/\tau}\right] \;=\; \frac{I(t)}{\tau}\, e^{t/\tau},

and since " $t$ " was just a placeholder for “any time,” we are free to relabel it. Call the dummy variable $s$ — the equation is identically true for $s$ in place of $t$ :

\frac{d}{ds}\!\left[V(s)\, e^{s/\tau}\right] \;=\; \frac{I(s)}{\tau}\, e^{s/\tau}.

Now $t$ stays as the fixed lower endpoint and $t + \Delta t$ stays as the fixed upper endpoint, while $s$ runs across the interval $[t,\, t + \Delta t]$ . The integral $\int_t^{t+\Delta t} \cdots\, ds$ is unambiguous, and the answer it produces will be a function of the endpoints $t$ and $t + \Delta t$ — which is exactly what we want, since the goal is to compute $V(t + \Delta t)$ from $V(t)$ .

Aside — the fundamental theorem of calculus

The fundamental theorem of calculus says that if $F$ is any antiderivative of $f$ (i.e. $F'(s) = f(s)$ ), then

\int_a^b f(s)\, ds \;=\; F(b) - F(a).

In words: integrating a derivative recovers the difference of the endpoint values. We don’t need to know how the integrand behaves anywhere except at the two endpoints — the antiderivative already captures everything in between. This is exactly the situation we have: the LHS of (10a) is itself a derivative of $V(s)\,e^{s/\tau}$ , so its antiderivative is just $V(s)\,e^{s/\tau}$ itself, and the integral collapses to the difference between its values at $s = t + \Delta t$ and $s = t$ :

\int_t^{t+\Delta t} \frac{d}{ds}\!\left[V(s)\, e^{s/\tau}\right] ds = V(t + \Delta t)\, e^{(t + \Delta t)/\tau} - V(t)\, e^{t/\tau}.

Equating to the integrated right-hand side:

V(t + \Delta t)\, e^{(t + \Delta t)/\tau} - V(t)\, e^{t/\tau} = \int_t^{t+\Delta t} \frac{I(s)}{\tau}\, e^{s/\tau}\, ds \tag{10b}

Step 3 — isolate $V(t + \Delta t)$

Divide both sides of (10b) by $e^{(t+\Delta t)/\tau}$ :

V(t + \Delta t) - V(t)\, e^{t/\tau}\, e^{-(t+\Delta t)/\tau} = \frac{1}{\tau}\, e^{-(t+\Delta t)/\tau} \int_t^{t+\Delta t} I(s)\, e^{s/\tau}\, ds

The decay coefficient on $V(t)$ simplifies using the exponent rule $e^{a}\, e^{b} = e^{a+b}$ :

e^{t/\tau}\, e^{-(t+\Delta t)/\tau} \;=\; e^{(t - t - \Delta t)/\tau} \;=\; e^{-\Delta t/\tau}.

That gives the homogeneous part of the solution: $V(t + \Delta t) = e^{-\Delta t/\tau}\, V(t) + \text{(integral term)}$ .

For the integral term, the factor $e^{-(t+\Delta t)/\tau}$ sits outside the integral but is constant with respect to the integration variable $s$ — neither $t$ nor $\Delta t$ depend on $s$ , they are fixed bounds. A standard property of integrals is that any factor independent of the integration variable can be moved freely across the integral sign:

\alpha \int_t^{t+\Delta t} f(s)\, ds \;=\; \int_t^{t+\Delta t} \alpha\, f(s)\, ds, \qquad \alpha \text{ independent of } s.

Applying that here with $\alpha = e^{-(t+\Delta t)/\tau}$ and $f(s) = I(s)\, e^{s/\tau}$ :

e^{-(t+\Delta t)/\tau} \int_t^{t+\Delta t} I(s)\, e^{s/\tau}\, ds \;=\; \int_t^{t+\Delta t} I(s)\, e^{s/\tau}\, e^{-(t+\Delta t)/\tau}\, ds.

Inside the integrand the two exponentials now multiply each other. Combining via the same exponent rule:

e^{s/\tau}\, e^{-(t+\Delta t)/\tau} \;=\; e^{(s - t - \Delta t)/\tau}.

So the integral term becomes $\frac{1}{\tau} \int_t^{t+\Delta t} I(s)\, e^{(s - t - \Delta t)/\tau}\, ds$ , and we can write the full update as

V(t + \Delta t) = e^{-\Delta t/\tau}\, V(t) + \frac{1}{\tau} \int_t^{t+\Delta t} I(s)\, e^{(s - t - \Delta t)/\tau}\, ds

Step 4 — change of variable $\xi = s - t$

Then $ds = d\xi$ , $s = t + \xi$ , and $s - t - \Delta t = \xi - \Delta t$ . The limits transform as $s = t \Rightarrow \xi = 0$ and $s = t + \Delta t \Rightarrow \xi = \Delta t$ .

Aside — why

ds = d\xi

The substitution defines a new variable $\xi$ as a function of $s$ via $\xi = s - t$ . Inside the integral, $s$ ranges over $[t,\, t + \Delta t]$ and $t$ is the fixed lower endpoint (recall the aside on renaming) — $t$ does not vary as $s$ moves through the interval. So $\xi$ is a function of just $s$ , with derivative

\frac{d\xi}{ds} \;=\; \frac{d}{ds}\,(s - t) \;=\; 1 - 0 \;=\; 1.

The differentials are then related by $d\xi = (d\xi/ds)\, ds = 1 \cdot ds = ds$ . Equivalently: an infinitesimal step in $s$ produces an equal infinitesimal step in $\xi$ , because the substitution is a pure shift (translation by $-t$ ) with slope 1.

This is what licenses the $u$ -substitution rule $\int f(\xi(s))\,(d\xi/ds)\, ds = \int f(\xi)\, d\xi$ ; with $d\xi/ds = 1$ the Jacobian factor is trivial and we can just relabel the variable directly.

Substituting:

V(t + \Delta t) = e^{-\Delta t/\tau}\, V(t) + \frac{1}{\tau} \int_0^{\Delta t} I(t + \xi)\, e^{-(\Delta t - \xi)/\tau}\, d\xi \tag{10}

This is exact for any forcing $I(t)$ — no approximation has been made yet.

Step 5 — replace the Dirac by a rectangle

So far we have an exact equation (10) for $V_{t+1}$ given $V_t$ and the forcing $I$ over the step. Now we hit the problem the rectangle approximation exists to solve. Suppose exactly one spike lands inside the step at some time $t_k \in [t,\, t + \Delta t]$ . The Dirac forcing for that spike alone is

I(t + \xi) \;=\; W\, \delta\!\bigl(t + \xi - t_k\bigr),

i.e. the forcing is zero everywhere across the step except at $\xi = t_k - t$ , where it is an infinitely tall pulse of area $W$ . Plugging that into the integral in (10),

\frac{1}{\tau} \int_0^{\Delta t} e^{-(\Delta t - \xi)/\tau}\, W\, \delta\!\bigl(t + \xi - t_k\bigr)\, d\xi,

and using the sifting property of the Dirac — $\int f(\xi)\, \delta(\xi - a)\, d\xi = f(a)$ whenever $a$ lies inside the integration interval — the integrand is collapsed to its value at the spike’s location, $\xi = t_k - t$ . The contribution to $V_{t+1}$ from that single spike comes out as

\frac{W}{\tau}\, e^{-(\Delta t - (t_k - t))/\tau}.

Read the exponent: $\Delta t - (t_k - t)$ is the amount of time between the spike and the end of the step. So the spike contributes more if it lands late in the step (small remaining time means little decay before reaching $t + \Delta t$ ) and less if it lands early (long remaining time means heavy exponential decay). In other words: the exact answer depends on where inside the step the spike arrived — on the sub-step arrival time $t_k - t$ , a continuous number between $0$ and $\Delta t$ .

That sub-step time is information the discrete simulator does not carry. Once we have agreed to integrate in steps of size $\Delta t$ , the only data we keep per neuron per step is “did at least one spike arrive in this step?” — a single Boolean $s_t \in \{0, 1\}$ . We have collapsed the continuous spike time onto a binary step indicator, and the moment we do that we forfeit the dependence on $t_k - t$ that the exact integral wants.

Aside — why the bias

b

doesn't appear in that contribution

The full forcing from (9) is $I(t) = \sum_k W\, \delta(t - t_k) + b$ — a Dirac train plus a constant bias. The integral in (10) is linear in $I$ , so the two parts contribute independently and we can write

\int_0^{\Delta t} e^{-(\Delta t - \xi)/\tau}\, I(t + \xi)\, d\xi \;=\; \underbrace{\int_0^{\Delta t} e^{-(\Delta t - \xi)/\tau}\, W\, \delta(t + \xi - t_k)\, d\xi}_{\text{spike part}} \;+\; \underbrace{\int_0^{\Delta t} e^{-(\Delta t - \xi)/\tau}\, b\, d\xi}_{\text{bias part}}.

The expression $(W/\tau)\, e^{-(\Delta t - (t_k - t))/\tau}$ shown above is just the spike part, isolated because that is the term whose value depends on the unknown sub-step arrival time $t_k - t$ . The bias part is a clean integral of a constant against the same kernel and evaluates to $b\,(1 - \beta)$ regardless of where the spike falls (it doesn’t depend on $t_k$ at all). We keep the spike part separate here only because that is what needs the rectangle approximation; the bias rides along untouched and re-enters in (11) below.

We replace each Dirac that lands in $[t, t+\Delta t]$ by a rectangle of equal area $W$ and width $\Delta t$ , i.e. constant height $W/\Delta t$ . Writing $s_t \in \{0, 1\}$ for the indicator that a spike occurred in step $t$ , the effective forcing over the step is

I(t + \xi) \;\to\; \frac{W\, s_t}{\Delta t} + b, \qquad \xi \in [0, \Delta t]

The bias $b$ is already constant, so it is unmodified. Define $\beta \equiv e^{-\Delta t / \tau}$ and $V_t \equiv V(t)$ . Substituting into (10):

V_{t+1} = \beta\, V_t + \frac{1}{\tau} \int_0^{\Delta t} e^{-(\Delta t - \xi)/\tau} \left[\frac{W\, s_t}{\Delta t} + b\right] d\xi \tag{11}

Step 6 — pull the constants out of the integral

The bracket is independent of $\xi$ , so

V_{t+1} = \beta\, V_t + \left[\frac{W\, s_t}{\Delta t} + b\right] \cdot \underbrace{\frac{1}{\tau} \int_0^{\Delta t} e^{-(\Delta t - \xi)/\tau}\, d\xi}_{\text{call this } J}

Step 7 — evaluate $J$ by $u$ -substitution

Let $u = \Delta t - \xi$ . Then $du = -d\xi$ , i.e. $d\xi = -du$ . The limits transform as $\xi = 0 \Rightarrow u = \Delta t$ and $\xi = \Delta t \Rightarrow u = 0$ . Substituting:

J = \frac{1}{\tau} \int_{u=\Delta t}^{u=0} e^{-u/\tau}\, (-du)

The negative sign flips the limits:

J = \frac{1}{\tau} \int_0^{\Delta t} e^{-u/\tau}\, du

The antiderivative of $e^{-u/\tau}$ with respect to $u$ is $-\tau\, e^{-u/\tau}$ (check: $\frac{d}{du}[-\tau\, e^{-u/\tau}] = -\tau \cdot (-1/\tau)\, e^{-u/\tau} = e^{-u/\tau}$ ). Evaluate:

J = \frac{1}{\tau}\, \Big[ -\tau\, e^{-u/\tau} \Big]_0^{\Delta t} = \frac{1}{\tau}\, \Big[ -\tau\, e^{-\Delta t/\tau} - (-\tau \cdot 1) \Big] = \frac{1}{\tau}\, \big[\, \tau - \tau\, e^{-\Delta t/\tau} \,\big] = 1 - e^{-\Delta t/\tau} = 1 - \beta \tag{12}

Step 8 — substitute $J$ back and distribute

From Step 6 with $J = 1 - \beta$ :

V_{t+1} = \beta\, V_t + (1 - \beta) \left[\frac{W\, s_t}{\Delta t} + b\right] = \beta\, V_t + \frac{1 - \beta}{\Delta t}\, W\, s_t + (1 - \beta)\, b

Renaming $V \to U$ to match the code’s symbol for membrane potential:

\boxed{U_{t+1} = \beta\, U_t + \frac{1 - \beta}{\Delta t}\, W\, s_t + (1 - \beta)\, b} \tag{13}

In the code, the two prefactors are computed once per step as spike_scale $= (1-\beta)/\Delta t$ and bias_scale $= 1 - \beta$ ; see models.py:705-710.

Step 9 — $\Delta t$ -invariance, part 1: bias steady state

In the absence of spikes ( $s_t = 0$ ), (13) reduces to $U_{t+1} = \beta\, U_t + (1-\beta)\, b$ . Setting $U_{t+1} = U_t = \bar U$ for the fixed point:

\bar U = \beta\, \bar U + (1 - \beta)\, b \;\;\Longrightarrow\;\; \bar U - \beta\, \bar U = (1 - \beta)\, b \;\;\Longrightarrow\;\; (1 - \beta)\, \bar U = (1 - \beta)\, b \;\;\Longrightarrow\;\; \bar U = b

The steady-state membrane equals $b$ regardless of $\Delta t$ . ✓

Step 10 — $\Delta t$ -invariance, part 2: single-spike time integral

Consider a network at rest, $U_0 = 0$ , that receives one spike at step 0 ( $s_0 = 1$ , $s_t = 0$ for $t \ge 1$ ), with $b = 0$ . Then from (13):

U_1 = \beta \cdot 0 + \frac{1 - \beta}{\Delta t}\, W \cdot 1 + 0 = K, \quad K \equiv \frac{1 - \beta}{\Delta t}\, W

For $t \ge 1$ there is no further forcing, so $U_{t+1} = \beta\, U_t$ , which means $U_t = K\, \beta^{t-1}$ for $t \ge 1$ . To get a $\Delta t$ -invariant quantity, integrate $U$ against time using a left-Riemann sum with step $\Delta t$ :

\int_0^\infty U(t)\, dt \;\approx\; \sum_{t=0}^\infty U_t\, \Delta t = U_0\, \Delta t + \sum_{t=1}^\infty K\, \beta^{t-1}\, \Delta t = 0 + K\, \Delta t \sum_{m=0}^\infty \beta^m

The sum is a geometric series with ratio $\beta \in (0, 1)$ , so $\sum_{m=0}^\infty \beta^m = 1/(1 - \beta)$ . Substituting:

\int_0^\infty U(t)\, dt \;\approx\; \frac{K\, \Delta t}{1 - \beta} = \frac{1}{1 - \beta} \cdot \frac{1 - \beta}{\Delta t}\, W \cdot \Delta t = W

The result is exactly $W$ , independent of $\Delta t$ . The corresponding continuous-time impulse response of (9) is $V(t) = (W/\tau)\, e^{-t/\tau}\, \mathbf{1}_{t \ge 0}$ , whose time integral is

\int_0^\infty \frac{W}{\tau}\, e^{-t/\tau}\, dt = \frac{W}{\tau} \cdot \tau \,\big[\,1 - 0\,\big] = W

The discrete and continuous time integrals agree. ✓

Step 11 — small- $\Delta t$ limit (sanity check)

Taylor-expand $\beta = e^{-\Delta t/\tau}$ around $\Delta t = 0$ :

\beta = 1 - \frac{\Delta t}{\tau} + \frac{1}{2}\!\left(\frac{\Delta t}{\tau}\right)^{\!2} - \cdots \;\;\Longrightarrow\;\; 1 - \beta = \frac{\Delta t}{\tau} - \frac{1}{2}\!\left(\frac{\Delta t}{\tau}\right)^{\!2} + \cdots

To leading order in $\Delta t/\tau$ :

$(1 - \beta)/\Delta t = 1/\tau - \Delta t/(2\tau^2) + O(\Delta t^2) \to 1/\tau$ .
$(1 - \beta) = \Delta t/\tau + O(\Delta t^2)$ , so $(1 - \beta)/\Delta t = 1/\tau + O(\Delta t)$ .

In words: the per-spike kick $(1-\beta)/\Delta t \cdot W$ approaches $W/\tau$ — i.e. the impulse response of the continuous LIF — and the per-step bias contribution $(1-\beta)\, b$ shrinks linearly with $\Delta t$ . The number of steps per ms scales as $1/\Delta t$ , so the per-millisecond bias contribution is $(1-\beta)\, b / \Delta t \to b/\tau$ , also $\Delta t$ -invariant.

Contrast with snnTorch

In the snnTorch discretisation, $W$ and $b$ have prefactor $1$ regardless of $\Delta t$ , so the per-spike kick is $W$ (not $W/\tau$ ) and the per-step bias is $b$ (not $b\, \Delta t/\tau$ ). Changing $\Delta t$ at eval time therefore shifts the effective input scale in standard-snn, while ZOH leaves it fixed — verified empirically in 013 — dt sensitivity.

Reset modes

Two choices on spike, controlled by —reset-mode:

zero — hard reset $U_{t+1} \leftarrow 0$ . Discards overshoot. Standard for biological neurons.
subtract — $U_{t+1} \leftarrow U_{t+1} - \theta$ . Preserves overshoot above threshold for the next step. Standard for snnTorch tutorials.

Defaults: standard-snn uses subtract under tutorial-mode; cuba uses zero.

Refractory period

Optional, controlled by —ref-ms. Implemented as a per-neuron can_fire mask in snn_lif_step: during refractory the membrane is clamped to $V_{\text{reset}}$ and the spike output is forced to zero. Default 0 ms (no refractory) for the CUBA-family rungs.

Where this leads

The two CUBA rungs span the design space of a single LIF cell with current-based inputs. The next step up the ladder adds genuine conductance-based synapses and a recurrent E-I loop — the COBA and PING models, documented in COBANet and exercised throughout the notebooks.