Clock driven synapse

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! :blush:
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. :blush:

1 Like