fix invalid refractory period of reduced models with Delta synapses

[NeoNeuron]

Description

How Has This Been Tested

The setting of refractory period in reduced models is invalid when I using Delta synapses. Here is my test code.

import brainpy.math as bm
import brainpy as bp
import matplotlib.pyplot as plt
bm.set_platform('cpu')

class delta_net(bp.Network):
  def __init__(self, method='exp_auto', **kwargs):
    super(delta_net, self).__init__(**kwargs)

    # neurons
    pars = dict(V_rest=-60., V_th=-50., V_reset=-60., tau=20., tau_ref=20.)
    E = bp.neurons.LIF(1, **pars, method=method)

    # synapses
    w_e = 1.2  # excitatory synaptic weight
    
    # ffwd Poisson drive
    self.Poisson = bp.neurons.PoissonGroup(1, freqs=1000, seed=0)
    self.P2E = bp.synapses.Delta(self.Poisson, E, bp.conn.One2One(), g_max=w_e)

    self.E = E
    
net = delta_net()

T = 100
runner = bp.DSRunner(
    net, monitors=['E.spike', 'E.V'], dt=0.1)
t = runner.run(T)
print('mean firing rate:', runner.mon['E.spike'].sum()/T*1e3, 'Hz')

# visualization
fig, ax = plt.subplots(2, 1, figsize=(6, 4), sharex='all')
bp.visualize.raster_plot(runner.mon.ts, runner.mon['E.spike'],
                         title='Spikes of Excitatory Neurons', 
                         marker='|', markersize=20, color='r', ax=ax[0])
ax[0].set_xlim(0,T)
ax[0].set_ylim(-.5,.5)

ax[1].plot(runner.mon.ts, runner.mon['E.V'][:,0])
ax[1].set_ylabel('Membrane Potential (mV)')
ax[1].axhline(-50, ls='--', c='r')
ax[1].axhline(-60, ls='--', c='r')
fig.savefig('delta_net.png', dpi=300, bbox_inches='tight')

Typically since the tau_ref=20., the maximum firing rate of the LIF neuron should be 50 Hz. However, the simulation above returns the 110 Hz firing rate.

The reason of this issue is that Delta synapses directly target to the membrane potentials self.V.value, which are updated before the integration of membrane potential.

BrainPy/brainpy/_src/neurons/reduced_models.py

Line 431 in d2c40da

V = self.integral(self.V.value, t, x, dt)

Therefore, the neuron will integrate the driven from Delta synapses even thought it's in refractory period.

https://github.com/brainpy/BrainPy/blob/d2c40daabbef8acbe82cc93e31e5bee42735c855/brainpy/_src/neurons/reduced_models.py#LL438C40-L438C40

I've change the update of V according to refractory with V = bm.where(refractory, self.V_reset, V). It solves the issue of invalidness of refractory.

However, I'm not sure whether my change will cause any trouble for networks in training mode.

Types of changes

• Bug fix (non-breaking change which fixes an issue)

Checklist

• Code follows the code style of this project.
• Changes follow the CONTRIBUTING guidelines.
• Update necessary documentation accordingly.
• Lint and tests pass locally with the changes.
• Check issues and pull requests first. You don't want to duplicate effort.

Other information

Two typos in the tutorial documents are fixed.

[chaoming0625]

It's a great PR. We will discuss more for solving this issue.

[chaoming0625]

Actually, running the model with brainpy.synapse.Delta needs the following setting:

class delta_net(bp.Network):
  def __init__(self, method='exp_auto', **kwargs):
    super(delta_net, self).__init__(**kwargs)

    # neurons
    pars = dict(V_rest=-60., V_th=-50., V_reset=-60., tau=20., tau_ref=20.)
    E = bp.neurons.LIF(1, **pars, method=method, ref_var=True)

    # synapses
    w_e = 1.2  # excitatory synaptic weight

    # ffwd Poisson drive
    self.Poisson = bp.neurons.PoissonGroup(1, freqs=1000, seed=0)
    self.P2E = bp.synapses.Delta(self.Poisson,
                                 E,
                                 bp.conn.One2One(), g_max=w_e,
                                 post_ref_key='refractory')

    self.E = E

The key difference is:

1. initializing the neuron model with brainpy.neuron.LIF(..., ref_var=True).
2. initializing the synapse model with bp.synapses.Delta(..., post_ref_key='refractory').

With these setting, models with refractory periods will produce the right results.

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

bm.set_platform('cpu')


class delta_net(bp.Network):
  def __init__(self, method='exp_auto', **kwargs):
    super(delta_net, self).__init__(**kwargs)

    # neurons
    pars = dict(V_rest=-60., V_th=-50., V_reset=-60., tau=20., tau_ref=20.)
    E = bp.neurons.LIF(1, **pars, method=method, ref_var=True)

    # synapses
    w_e = 1.2  # excitatory synaptic weight

    # ffwd Poisson drive
    self.Poisson = bp.neurons.PoissonGroup(1, freqs=1000, seed=0)
    self.P2E = bp.synapses.Delta(self.Poisson,
                                 E,
                                 bp.conn.One2One(), g_max=w_e,
                                 post_ref_key='refractory')

    self.E = E


net = delta_net()

T = 100
runner = bp.DSRunner(
  net, monitors=['E.spike', 'E.V'], dt=0.1)
t = runner.run(T)
print('mean firing rate:', runner.mon['E.spike'].sum() / T * 1e3, 'Hz')

# visualization
fig, ax = plt.subplots(2, 1, figsize=(6, 4), sharex='all')
bp.visualize.raster_plot(runner.mon.ts, runner.mon['E.spike'],
                         title='Spikes of Excitatory Neurons',
                         marker='|', markersize=20, color='r', ax=ax[0])
ax[0].set_xlim(0, T)
ax[0].set_ylim(-.5, .5)

ax[1].plot(runner.mon.ts, runner.mon['E.V'][:, 0])
ax[1].set_ylabel('Membrane Potential (mV)')
ax[1].axhline(-50, ls='--', c='r')
ax[1].axhline(-60, ls='--', c='r')
# fig.savefig('delta_net.png', dpi=300, bbox_inches='tight')
plt.show()

[chaoming0625]

Sorry, I have posted the wrong image of the running. The following figure is the correct one:

[NeoNeuron]

Thanks for reply. The settings you presented fixed my issue well. Sorry that I just miss the post_ref_key part in the Delta synapses. However, I think it's less apparent to notice this special setting for Delta synapses and neurons with refractory period. And since the I believe Delta synapse is widely used in the simulation of large balanced network, maybe we can add a tutorial for this or a few more words in the docs of brainpy.synapses.Delta object.

[ztqakita]

Thank you for a very good suggestion, we will add clearer notes as soon as possible to help people understand the model better.

[ztqakita]

I believe this issue has been well addressed and further updates will be coming soon.