Description of problem
I am trying to implement the homeostatic inhibitory STDP rule used in Zenke (2015) in Brian2. The rule uses a symmetric STDP kernel that is modulated in size and sign by a global factor G, which corresponds to the difference between a (smoothed) estimate of the population firing rate H and a desired firing rate gamma. I am trying to implement this rule using a synapse model with a linked variable. However, my current attempt results in a cryptic error message.
Full traceback of error (if relevant)
Traceback (most recent call last):
File "/home/paul/work/gothard/code/inhibitory_homeostatic_plasticity_with_shared_global_factor.py", line 61, in <module>
synapses = b2.Synapses(
^^^^^^^^^^^^
File "/home/paul/opt/miniconda3/envs/brian2/lib/python3.11/site-packages/brian2/synapses/synapses.py", line 853, in __init__
model.check_flags(
File "/home/paul/opt/miniconda3/envs/brian2/lib/python3.11/site-packages/brian2/equations/equations.py", line 1234, in check_flags
raise ValueError(
ValueError: Equations of type 'parameter' cannot have a flag 'linked', only the following flags are allowed: ['constant', 'shared']
Minimal code to reproduce problem
"""
This script implements a STDP rule with a global/shared factor as in Zenke et al. (2015).
"""
import numpy as np
import matplotlib.pyplot as plt
import brian2 as b2
b2.start_scope()
net = b2.Network(b2.collect())
# set (target) population firing rate;
# in a non-toy simulation, spikes would be recorded from a large pyramidal neuron population
total_epochs = 11
epoch_length = 1 * b2.second
population_size = 1
spike_rates = np.geomspace(1, 50, total_epochs)
spike_rates_array = b2.TimedArray(spike_rates * b2.Hz, dt=epoch_length)
spike_generators = b2.PoissonGroup(population_size, rates="spike_rates_array(t)")
net.add(spike_generators)
# compute global factor G
gamma = 5 * b2.Hz
tau_H = 10 * b2.second
global_factor_model = """
G = H - gamma : Hz
dH/dt = -H/tau_H : Hz
"""
global_factor = b2.NeuronGroup(1, model=global_factor_model)
net.add(global_factor)
spike_aggregator = b2.Synapses(spike_generators, global_factor, on_pre="H_post += 1/(tau_H * population_size)")
spike_aggregator.connect()
net.add(spike_aggregator)
# initialize the pre- and post-synaptic neuron
indices = np.zeros(total_epochs)
times = epoch_length * np.arange(total_epochs) + 0.5 * epoch_length
presynaptic_neuron = b2.SpikeGeneratorGroup(1, indices, times)
net.add(presynaptic_neuron)
delta = 10 * b2.msecond
postsynaptic_neuron = b2.SpikeGeneratorGroup(1, indices, times + delta)
net.add(postsynaptic_neuron)
# initialize the synapse
tau_stdp = 20. * b2.msecond
eta = 1. * b2.nsiemens * b2.second
synapse_model = """
G : Hz (linked)
wij : siemens
dzi/dt = -zi / tau_stdp : 1 (event-driven)
dzj/dt = -zj / tau_stdp : 1 (event-driven)
"""
on_pre = """
zi += 1.
wij += eta * G * zj
"""
on_post = """
zj += 1.
wij += eta * G * (zi + 1)
"""
synapses = b2.Synapses(
presynaptic_neuron, postsynaptic_neuron,
model=synapse_model, on_pre=on_pre, on_post=on_post
)
synapses.connect(i=0, j=0)
# synapses.G = b2.linked_var(global_factor.G)
synapses.G = b2.linked_var(global_factor, "G")
synapses.wij[:] = 0 * b2.nsiemens
net.add(synapses)
# monitor synaptic variables
monitor = b2.StateMonitor(synapses, ["G", "wij"], record=True, dt=epoch_length)
net.add(monitor)
# run
net.run(total_epochs * epoch_length)
# plot
fig, axes = plt.subplots(2, 1, sharex=True)
axes[0].plot(np.squeeze(monitor.G / b2.Hz), color="gray")
axes[0].set_ylabel("Global factor [Hz]")
axes[1].plot(np.diff(np.squeeze(monitor.wij / b2.nsiemens)), color="gray")
axes[1].set_ylabel("Weight [nS]")
axes[-1].set_xlabel("Epoch")
plt.show()
Any insight or workaround would be appreciated.
). The reason for the error is that the