Is it possible to implement a clock-driven synapse in Brian2?
I want to model a synapse where each neuron in a NeuronGroup has a variable ( r ) that evolves according to:
\frac{dr}{dt} = -\frac{r}{\tau} + u \delta(t)
where ( \delta(t) ) represents the spike train of the neuron.
Additionally:
- When a neuron spikes (i.e., crosses its threshold), ( r ) increases by ( u ).
- The synaptic input to the postsynaptic neuron is proportional to the presynaptic neuron’s ( r ), multiplied by the synaptic weight.
Would this be feasible in Brian2 to have such synapse? I am fine with defining the variable r but i dont know how i can update membrane potential at each time step based on r values and corresponding synaptic weight.
Hi @elnaz91. I am a bit confused, you are saying that \delta(t) represents “the spike train of the neuron” and that when a"a neuron spikes (i.e., crosses its threshold), ( r ) increases by ( u ).". But do you really mean that a neuron’s r should increase when its own r goes over the threshold?
Regarding the second part, if you want a synapse to provide continuous input to a post-synaptic neuron based on r of a pre-synaptic neuron, you can use a “summed variable” (see Synapses — Brian 2 2.8.0.4 documentation) – it sums over the given expression for each synapse connecting to the post-synaptic cell. So if your post-synaptic cell has a variable syn_inp
to denote the total synaptic input at every time step, you’ll write the following equation for the Synapses
:
syn_eqs = '''
w : 1 (constant)
syn_inp_post = r_pre * w : 1 (summed)
'''
Hope that helps!
1 Like
hat is correct. Since r represents the firing rates of each neuron and has a slower time constant compared to the membrane potential, I believe you have answered my question. Thanks a million!
Below is a simpler version of my code.
from brian2 import *
# Set the simulation parameters
start_scope()
tau = 20 * ms
tau_r = 100 * ms
tau_fast = 1 * ms
El = -80 * mV # Resting potential
Vr = -80 * mV # Reset potential
Vt = -30 * mV # Threshold potential
# Neuron model equations
eqs_neurons = '''
dv/dt = (El - v) / tau + I_fast + I_slow : volt
dr/dt = -r / tau_r : 1
dI_fast/dt = -I_fast / tau_fast : volt/second # Fast synaptic current decays exponentially
I_slow : volt/second # Slow synaptic current
'''
# Number of neurons
N = 10
M = 100 # Number of Poisson input sources
# Neuron group
neurons = NeuronGroup(N, eqs_neurons, threshold='v > Vt', reset='v = Vr; r += 1/tau_r*ms', refractory=2*ms, method='euler')
neurons.v = El # Initialize membrane potential
neurons.r = 0 # Initialize variable r
neurons.I_fast = 0
neurons.I_slow = 0
# Synapse equations for slow synapses
eqs_synapses = '''
w : 1 # Synaptic weight (dimensionless)
I_slow_post = w * r_pre * mV/ms : volt/second (summed)
'''
# Create synapses with random weights between 0 and 2
synapses = Synapses(neurons, neurons, model=eqs_synapses, method='euler')
synapses.connect()
synapses.w = rand()*3 # Random weights between 0 and 2 for stability
# Poisson input synapses for fast currents
P = PoissonGroup(M, rates=80*Hz)
input_synapses = Synapses(P, neurons, model='w : 1 (constant)', on_pre='I_fast_post += w *mV/ms')
input_synapses.connect()
input_synapses.w = rand()
# Monitors
state_monitor = StateMonitor(neurons, ('v', 'r', 'I_slow', 'I_fast'), record=True)
spike_monitor = SpikeMonitor(neurons)
# Run the simulation
run(2*second)
# Plotting
figure(figsize=(20, 12))
# Plot membrane potential of all neurons
subplot(221)
for i in range(N):
plot(state_monitor.t/ms, state_monitor.v[i]/mV, label=f'Neuron {i}')
xlabel('Time (ms)')
ylabel('Membrane Potential (mV)')
title('Membrane Potential of All Neurons')
legend(loc='upper right')
# Plot r variable for all neurons
subplot(222)
for i in range(N):
plot(state_monitor.t/ms, state_monitor.r[i], label=f'Neuron {i}')
xlabel('Time (ms)')
ylabel('r')
title('Variable r of All Neurons')
legend(loc='upper right')
# Plot I_fast and I_slow currents for all neurons
subplot(223)
for i in range(1):
plot(state_monitor.t/ms, state_monitor.I_fast[i]*ms/mV, label=f'I_fast Neuron {i}')
xlabel('Time (ms)')
ylabel('I_fast (mV)')
title('Fast Synaptic Currents')
legend(loc='upper right')
subplot(224)
for i in range(1):
plot(state_monitor.t/ms, state_monitor.I_slow[i]*ms/mV, label=f'I_slow Neuron {i}')
xlabel('Time (ms)')
ylabel('I_slow (mV)')
title('Slow Synaptic Currents')
legend(loc='upper right')
# Show all current plots
show()
# Spike raster plot
figure(figsize=(10, 6))
plot(spike_monitor.t/ms, spike_monitor.i, 'k.', markersize=5)
xlabel('Time (ms)')
ylabel('Neuron index')
title('Spike Raster Plot')
show()
1 Like
Hi @elnaz91. Looks good to me! I misunderstood your model earlier, I thought that r would have been the only variable of your model (i.e. also represent the “membrane potential”) – but in your formulation it makes perfect sense.
1 Like
Thank you, Marcel! Your quick response means a lot to the Brian2 community.
1 Like