Difference between the recorded conductance of Alpha synapse and the convolution of alpha function to spike trains.

[CloudyDory]

The response of an Alpha synapse after a single spike input is:

g_syn(t) = g_max * t/tau * exp(-t/tau)

For a train of input spikes (represented by binary vector s(t) = [0 0 0 1 0 1 0 0 ...]), the response will be the convolution between s(t) and g_syn(t).

In BrainPy, for computational efficiency, the convolution is replaced by the following differential equations:

g_syn(t) = g_max * g
dg/dt = -g/tau + h
dh/dt = -h/tau
h <-- h + s(t)

However, the recorded g_syn(t) by this method differs by the convolution method by a factor of tau. The following code illustrates this:

import numpy as np
import brainpy as bp
import brainpy.math as bm
import matplotlib.pyplot as plt

dt = 0.1     # ms
bm.set_dt(dt)

neu1 = bp.neurons.LIF(1, method='rk4')
neu2 = bp.neurons.LIF(1, method='rk4')
syn1 = bp.synapses.Alpha(neu1, neu2, bp.connect.All2All(), g_max=1.0, tau_decay=8.6, output=bp.synouts.CUBA(), method='rk4')
net = bp.Network(pre=neu1, syn=syn1, post=neu2)

runner = bp.DSRunner(net, inputs=[('pre.input', 25.)], monitors=['pre.V', 'pre.spike', 'post.V', 'syn.g', 'syn.h'], dt=dt, jit=True)
runner.run(150.)

time = np.maximum(np.arange(0, 200, dt) - dt, 0)  # Brainpy lags one time point
kernel = syn1.g_max * time/syn1.tau_decay * np.exp(-time/syn1.tau_decay)
g_conv = np.convolve(runner.mon['pre.spike'].squeeze(), kernel, mode='full')[:runner.mon['pre.spike'].shape[0]] * syn1.tau_decay

plt.figure()
plt.plot(runner.mon.ts, runner.mon['syn.g'], label='g_recorded')
plt.plot(runner.mon.ts, g_conv, label='g_conv*tau_decay')
plt.legend()
plt.xlabel('Time (ms)')
plt.ylabel('Conductance')
plt.show()

The key is the h <-- h + s(t) part. If I change it to h <-- h + s(t)/tau, the outputs from differential equation form and convolution form will be the same. However, I have checked other simulator (https://brian2.readthedocs.io/en/stable/user/converting_from_integrated_form.html), but it seems that it uses similar implementation as in BrainPy (i.e. directly shifting h without divided by tau). So I am a bit confused. Is the larger response of differential equation form correct?

[chaoming0625]

Thanks for the report! This is a great issue!

Let me first guess the reason behind this error.

The error may result from the incorrect discretization of the Alpha synapse model.

Current discretization is

$$ \begin{aligned} &\frac{d g}{d t}=-\frac{g}{\tau_{\mathrm{decay}}}+h \\ &\frac{d h}{d t}=-\frac{h}{\tau_{\text {rise }}}+ \delta\left(t_{0}-t\right), \end{aligned} $$

However, it should be

$$ \begin{aligned} & \tau_{\mathrm{decay}} \frac{d g}{d t}=-g+h \\ &\frac{d h}{d t}=-\frac{h}{\tau_{\text {rise }}}+ \delta\left(t_{0}-t\right), \end{aligned} $$

I will recheck the math later.

[ztqakita]

Thanks for the report, we have rechecked the error and fixed the bugs in BrainPy. Please see #578 .