# Handling time-constant singularities in integrate-and-fire neurons¶

This notebook details how NEST handles the numerical instability of the exact integration propagator matrix $$P = e^{A h}$$ which arise if $$\tau_m\approx \tau_{\text{syn}}$$. For an overview over exact integration of integrate-and-fire neuron subthreshold dynamics, please see the exact integration documentation.

We illustrate the approach for neurons with alpha-shaped currents, where the synaptic current is described by two differential equations. For exponential-shaped currents, a similar but simpler treatment applies.

The singularity-handling code is implemented in libnestutil/iaf_propagator.[h,cpp].

## Preparations¶

We use SymPy to allow symbolic analysis of the propagator matrices and their limits.

[1]:

import sympy as sp
from sympy.matrices import zeros

sp.init_printing(use_latex=True)


Introduce formal variables for time constants, capacitance and time step $$h$$.

[2]:

tau_m, tau_s, C_m, h = sp.symbols("tau_m, tau_s, C_m, h", positive=True)


## The ODE matrix¶

The following matrix describes the ODE system for synaptic current and membrane potential (bottom row). It applies for singular and non-singular cases. For brevity, we write $$\tau_s$$ instead of $$\tau_{\text{syn}}$$.

[3]:

A = sp.Matrix([[-1 / tau_s, 0, 0], [1, -1 / tau_s, 0], [0, 1 / C_m, -1 / tau_m]])
A

[3]:

$\displaystyle \left[\begin{matrix}- \frac{1}{\tau_{s}} & 0 & 0\\1 & - \frac{1}{\tau_{s}} & 0\\0 & \frac{1}{C_{m}} & - \frac{1}{\tau_{m}}\end{matrix}\right]$

## Propagator in the non-singular case ($$\tau_m\neq \tau_s$$)¶

[4]:

P = sp.simplify(sp.exp(A * h))
P

[4]:

$\displaystyle \left[\begin{matrix}e^{- \frac{h}{\tau_{s}}} & 0 & 0\\h e^{- \frac{h}{\tau_{s}}} & e^{- \frac{h}{\tau_{s}}} & 0\\\frac{\tau_{m} \tau_{s} \left(- h \tau_{m} e^{\frac{h}{\tau_{m}}} + h \tau_{s} e^{\frac{h}{\tau_{m}}} - \tau_{m} \tau_{s} e^{\frac{h}{\tau_{m}}} + \tau_{m} \tau_{s} e^{\frac{h}{\tau_{s}}}\right) e^{- \frac{h}{\tau_{s}} - \frac{h}{\tau_{m}}}}{C_{m} \left(\tau_{m}^{2} - 2 \tau_{m} \tau_{s} + \tau_{s}^{2}\right)} & \frac{\tau_{m} \tau_{s} \left(- e^{\frac{h}{\tau_{m}}} + e^{\frac{h}{\tau_{s}}}\right) e^{- \frac{h \left(\tau_{m} + \tau_{s}\right)}{\tau_{m} \tau_{s}}}}{C_{m} \left(\tau_{m} - \tau_{s}\right)} & e^{- \frac{h}{\tau_{m}}}\end{matrix}\right]$
• The entries in the first two rows of $$P$$ are unproblematic, as is $$P_{33}$$.

• In the first two entries on the bottom row, $$P_{31}$$ and $$P_{32}$$, the denominators will vanish for $$\tau_s\to\tau_m$$.

• $$P_{32}$$ also appears in the propagator matrix for the case of exponential synaptic currents.

## Propagator in the singular case ($$\tau_s = \tau_m$$)¶

In this case, we have

[5]:

A_s = sp.Matrix([[-1 / tau_m, 0, 0], [1, -1 / tau_m, 0], [0, 1 / C_m, -1 / tau_m]])
A_s

[5]:

$\displaystyle \left[\begin{matrix}- \frac{1}{\tau_{m}} & 0 & 0\\1 & - \frac{1}{\tau_{m}} & 0\\0 & \frac{1}{C_{m}} & - \frac{1}{\tau_{m}}\end{matrix}\right]$

and the propagator becomes

[6]:

P_s = sp.simplify(sp.exp(A_s * h))
P_s

[6]:

$\displaystyle \left[\begin{matrix}e^{- \frac{h}{\tau_{m}}} & 0 & 0\\h e^{- \frac{h}{\tau_{m}}} & e^{- \frac{h}{\tau_{m}}} & 0\\\frac{h^{2} e^{- \frac{h}{\tau_{m}}}}{2 C_{m}} & \frac{h e^{- \frac{h}{\tau_{m}}}}{C_{m}} & e^{- \frac{h}{\tau_{m}}}\end{matrix}\right]$
• This is well-formed and non-singular.

• The “unproblematic” matrix elements agree with the non-singular case.

## Numeric stability of propagator elements¶

We will now show that the matrix elements of the non-singular case converge to those in the general case, so that we have overall

$\lim_{\tau_s\to\tau_m} P = P_s\;.$

Using symbolic algebra we find for $$\lim_{\tau_s\to\tau_m} P_{31} = P_{s,31}$$:

[7]:

P_31 = P.row(2).col(0)[0]
P_31_l = sp.limit(P_31, tau_s, tau_m)
P_31_l

[7]:

$\displaystyle \frac{h^{2} e^{- \frac{h}{\tau_{m}}}}{2 C_{m}}$

Test for mathematical equality as recommended in SymPy

[8]:

sp.simplify(P_31_l - P_s.row(2).col(0)[0]) == 0

[8]:

True

[9]:

P_32 = P.row(2).col(1)[0]
P_32_l = sp.limit(P_32, tau_s, tau_m)
P_32_l

[9]:

$\displaystyle \frac{h e^{- \frac{h}{\tau_{m}}}}{C_{m}}$
[10]:

sp.simplify(P_32_l - P_s.row(2).col(1)[0]) == 0

[10]:

True


## Approximation in the vicinity of the singularity¶

Since the propagator elements converge to the solution for the singular case, we can approximate the matrix elements near the singularity by expanding around $$\tau_m$$. We obtain

[11]:

P_31.series(x=tau_s, x0=tau_m, n=2)

[11]:

$\displaystyle \frac{h^{2} e^{- \frac{h}{\tau_{m}}}}{2 C_{m}} + \frac{h^{3} \left(- \tau_{m} + \tau_{s}\right) e^{- \frac{h}{\tau_{m}}}}{3 C_{m} \tau_{m}^{2}} + O\left(\left(- \tau_{m} + \tau_{s}\right)^{2}; \tau_{s}\rightarrow \tau_{m}\right)$
[12]:

P_32.series(x=tau_s, x0=tau_m, n=2)

[12]:

$\displaystyle \frac{h e^{- \frac{h}{\tau_{m}}}}{C_{m}} + \frac{h^{2} \left(- \tau_{m} + \tau_{s}\right) e^{- \frac{h}{\tau_{m}}}}{2 C_{m} \tau_{m}^{2}} + O\left(\left(- \tau_{m} + \tau_{s}\right)^{2}; \tau_{s}\rightarrow \tau_{m}\right)$

We thus have

\begin{align} P_{31} &= P_{s, 31} + \frac{2h}{3\tau_m^2}P_{s, 31}(\tau_s-\tau_m) + \mathcal{O}((\tau_s-\tau_m)^2)\\ P_{32} &= P_{s, 32} + \frac{h}{2\tau_m^2}P_{s, 32}(\tau_s-\tau_m) + \mathcal{O}((\tau_s-\tau_m)^2) \end{align}

Focusing on $$P_{32}$$ and dropping the quadratic term, we obtain

$\left|\frac{P_{32}-P_{s, 32}}{P_{s, 32}}\right|\approx \left|\frac{h(\tau_s-\tau_m)}{2\tau_m^2}\right| \ll 1$

where the inequality follows because $$|\tau_s-\tau_m|\ll \tau_m$$ by definition in the near-singular case and $$h<\tau_m$$ for all practical purposes.

Any violation of this inequality indicates numerical instability in the computation of $$P_{32}$$.

The corresponding inequality for $$P_{31}$$ is

$\left|\frac{P_{31}-P_{s, 31}}{P_{s, 31}}\right|\approx \left|\frac{2h(\tau_s-\tau_m)}{3\tau_m^2}\right| \ll 1\;.$

### Propagators rewritten as expressed in C++ implementation¶

• Implementation uses expm1(x) function which returns $$e^x-1$$.

• Show that expressions for $$P_{31}$$ and $$P_{32}$$ in implemenation are equivalent to expressions above.

[13]:

beta = (tau_m * tau_s) / (tau_m - tau_s)
gamma = beta / C_m
beta, gamma

[13]:

$\displaystyle \left( \frac{\tau_{m} \tau_{s}}{\tau_{m} - \tau_{s}}, \ \frac{\tau_{m} \tau_{s}}{C_{m} \left(\tau_{m} - \tau_{s}\right)}\right)$
[14]:

P_32i = gamma * sp.exp(-h / tau_s) * (sp.exp(h / tau_s - h / tau_m) - 1)
P_32i

[14]:

$\displaystyle \frac{\tau_{m} \tau_{s} \left(e^{\frac{h}{\tau_{s}} - \frac{h}{\tau_{m}}} - 1\right) e^{- \frac{h}{\tau_{s}}}}{C_{m} \left(\tau_{m} - \tau_{s}\right)}$
[15]:

sp.simplify(P_32 - P_32i) == 0

[15]:

True

[16]:

P_31i = gamma * sp.exp(-h / tau_s) * (beta * (sp.exp(h / tau_s - h / tau_m) - 1) - h)
P_31i

[16]:

$\displaystyle \frac{\tau_{m} \tau_{s} \left(- h + \frac{\tau_{m} \tau_{s} \left(e^{\frac{h}{\tau_{s}} - \frac{h}{\tau_{m}}} - 1\right)}{\tau_{m} - \tau_{s}}\right) e^{- \frac{h}{\tau_{s}}}}{C_{m} \left(\tau_{m} - \tau_{s}\right)}$
[17]:

sp.simplify(P_31 - P_31i) == 0

[17]:

True


### Numerical convergence experiments¶

• Compute propagator elements as implemented in code

• Test convergence against singular value for $$\tau_s\to\tau_m$$

• Test for different time steps $$h$$

• Very small time steps $$\mathcal{O}(1 \mu{s})$$ can occur in neurons with precise spike times and are thus relevant

• We set $$C_m=1$$ since it is just a scaling factor

[18]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import nest

nest.set_verbosity("M_ERROR")


-- N E S T --
Copyright (C) 2004 The NEST Initiative

Version: stinebuu_propagator_class@79327cc82
Built: Nov 29 2022 17:36:40

This program is provided AS IS and comes with
NO WARRANTY. See the file LICENSE for details.

Problems or suggestions?
Visit https://www.nest-simulator.org

Type 'nest.help()' to find out more about NEST.


[19]:

def calc_p(tau_m, tau_s, h):
beta = tau_s * tau_m / (tau_m - tau_s)
inv_beta = (tau_m - tau_s) / (tau_s * tau_m)
gamma = beta

p31 = gamma * np.exp(-h / tau_s) * (beta * np.expm1(h * inv_beta) - h)
p31s = 0.5 * h**2 * np.exp(-h / tau_m)

p32 = gamma * np.exp(-h / tau_s) * np.expm1(h * inv_beta)
p32s = h * np.exp(-h / tau_m)

return p31, p31s, p32, p32s

[20]:

tau_m = 10.0
tau_s = tau_m + np.logspace(-12, 0, 26)
h = np.array([[1, 0.1, 0.01, 0.001]]).T
p31, p31s, p32, p32s = calc_p(tau_m, tau_s, h)

fig = plt.figure(figsize=(12, 4))
ax1.loglog(np.abs(tau_s - tau_m), np.abs(p31 - p31s).T, "o-")
ax1.set_ylabel("$P - P_{singular}$")
ax1.set_ylim([1e-20, 5e-2])
ax1.set_xlabel(r"$\tau_s - \tau_m$ [ms]")
ax1.set_title("$P_{31}$ difference from singular value")
ax2.loglog(np.abs(tau_s - tau_m), np.abs(p32 - p32s).T, "o-", label=[f"h = {hv:.3f} ms" for hv in h[:, 0]])
ax2.set_ylim([1e-20, 5e-2])
ax2.set_xlabel(r"$\tau_s - \tau_m$ [ms]")
ax2.set_title("$P_{32}$ difference from singular value")
plt.legend();

• $$P_{32}$$ shows perfect convergence towards the singular value up to the limits of numerical accuracy

• This holds for all step sizes $$h$$

• This is plausible because $$P_{32}$$ contains multiplications only: gamma * np.exp(-h/tau_s) * np.expm1(h*inv_beta)

• The only difference is handled internally in expm1() and inv_beta goes to 0 in the limit

• Thus no singularity handling is needed for :math:P_{32}.

• $$P_{31}$$ converges only up to a point, which depends on the size of the time step $$h$$

• Numerical instability occurs for smaller differences in time constants

• Instability occurs earlier for smaller time steps

• The instability arises from the difference beta * np.expm1(h*inv_beta) - h

• There seems to be no way to reformulate the propagator to avoid this

#### $$P_{31}$$ instability and membrane time constant¶

[21]:

h = np.array([[1, 0.1, 0.01, 0.001]]).T
tau_ms = [0.1, 1, 10, 100]
delta_tau = np.logspace(-12, 0, 26)

fig = plt.figure(figsize=(12, 3))
for ix, tau_m in enumerate(tau_ms):
tau_s = tau_m + delta_tau
p31, p31s, _, _ = calc_p(tau_m, tau_s, h)

ax = fig.add_subplot(1, len(tau_ms), ix + 1)
l = ax.loglog(np.abs(tau_s - tau_m), np.abs(p31 - p31s).T, "o-", label=[f"h = {hv:.3f} ms" for hv in h[:, 0]])
ax.set_prop_cycle(None)
for hv in h[:, 0]:
ax.loglog(1e-8 * tau_m**2 / hv, 1e-16, "^")
ax.set_ylim([1e-20, 5e-2])
ax.set_xlabel(r"$\tau_s - \tau_m$ [ms]")
ax.set_title(f"$\\tau_m = {tau_m}$ ms")
if ix == 0:
ax.set_ylabel("$P - P_{singular}$")
plt.legend(loc="upper left")
else:
ax.set_yticklabels([]);

• The point at which numerical instability occurs clearly depends on the membrane time constant

• As indicated by the markers, the breakdown point is located roughly at

$\tau_s - \tau_m < 10^{-8}\times\frac{\tau_m^2}{h}$
• We can thus use

$(\tau_s - \tau_m)h < 10^{-8}\times \tau_m^2$

or

$h < 10^{-8}\times\frac{\tau_m^2}{|\tau_s - \tau_m|}$

as criterium: use $$P_{s, 31}$$ if this condition is fulfilled.

• To ensure some margin of safety, we can set the limit at $$10^{-7}\times\tau_m^2$$.

## Algorithm for propagator computing¶

1. Precompute values that can be precomputed independent of $$h$$, including useful inverses.

2. Compute $$P_{32}$$, which will always be stable; replace by singular limit only if a numerically non-normal or non-positive result occurs.

3. If $$h$$ is below stability limit (see inequality above), use singular $$P_{s, 31}$$, otherwise use $$P_{31}$$.

## Exploration¶

We will now show that the stability criterion explained above leads to a reasonable behavior for $$\tau_s\rightarrow\tau_m$$

### Simulation¶

• Create one neuron for each value of delta_tau

• Drive neurons with a single spike

• Measure resulting membrane potential

[22]:

tau_m = 10.0
h = 0.1
delta_tau = np.hstack(([0.0], np.logspace(-10, -1, 10)))

nest.ResetKernel()
nest.resolution = h

neurons = nest.Create("iaf_psc_alpha", n=len(delta_tau), params={"tau_m": tau_m, "tau_syn_ex": tau_m + delta_tau})
spike_gen = nest.Create("spike_generator", params={"spike_times": [1.0]})
vm = nest.Create("voltmeter", params={"interval": h})

nest.Connect(spike_gen, neurons, syn_spec={"weight": 100.0})
nest.Connect(vm, neurons)

nest.Simulate(10 * tau_m)

v = pd.DataFrame.from_records(vm.events).set_index("times")

[23]:

v.groupby("senders").V_m.plot(alpha=0.5, figsize=(8, 3))
plt.legend(delta_tau)
plt.xlabel("Time [ms]")
plt.ylabel("Membrane voltage [mV]")
plt.title("Response to single input spike for different tau_m - tau_s");


### Maximum of membrane potential¶

[24]:

V_max = v.groupby("senders").V_m.max()
plt.figure(figsize=(10, 3))
plt.subplot(1, 3, 1)
plt.semilogx([1e-10, 1], [V_max.iloc[0], V_max.iloc[0]], "--")
plt.semilogx(delta_tau[1:], V_max.iloc[1:], "o-", alpha=0.6)
plt.legend(("Singular limit", "$tau_s > tau_m$"))
plt.xlabel(r"$\tau_s - \tau_m$")
plt.ylabel("$V_{max} [mV]$")

plt.subplot(1, 3, 3)
plt.loglog(delta_tau[1:], V_max.iloc[1:] - V_max.iloc[0], "o-")
plt.xlabel(r"$\tau_s - \tau_m$")
plt.ylabel("$V_{max} - V_{max}^{singular} [mV]$");


The maximum membrane potential converges smoothly against the singular limit, indicating that no numerical instabilities occur.