Hello,
I am new to neuron modelling and Brian2, and I am starting a project that requires modifying the membrane potential equation of a neuron group. I need to model the membrane potential change, V(t) , based on the convolution of an exponential function with the spike trains S(t) from the neurons themselves. The effect of this convolution is scaled by a constant weight w . v += W . h*S
Could you please provide guidance or examples on how to model this in Brian2?
Thank you!
Elnaz
Hi @elnaz91. If I understand correctly, you want to have something like the mechanism in the Spike-Response model? I’d say the most straightforward way to model this is along these lines, where the v_ap variable represents the convolution of the spike train with an exponential kernel:
eqs = '''
v_final = v + v_ap : 1
dv/dt = ... : 1 # Diff eq for the membrane potential
dv_ap/dt = -v_ap/tau_ap : 1 # exponential kernel
'''
group = NeuronGroup(..., eqs,
threshold="v_final > ...",
reset="v_ap -= w" # Jump down with each spike
)
Hope this helps!
Thank you, @mstimberg, for your response. That is really helpful.
I also need to define it in a way that could record excitatory and inhibitory currents to each neuron within the network. I have two types of excitatory currents (one is this one, and the other is the usual excitatory current that increases membrane potential by predefined values) and one is inhibitory. what do you think if I write membrane eqs as
eqs_neurons = ‘’’
dv/dt = (-v + I_exc + I_inh) / tau_membrane : volt
I_exc : volt
I_inh : volt
‘’’
eqs_exc1_synapse = ‘’’
w : volt
dx/dt = -x / tau_exc1 : 1 (event-driven)
I_exc_post += w * x : volt
‘’’
eqs_exc2_synapse = ‘’’
w : volt
‘’’
eqs_inh_synapse = ‘’’
w : volt
‘’’
exc1_syn = Synapses(neurons, neurons, model=eqs_exc1_synapse, on_pre=‘x += 1’)
exc1_syn.connect(p=0.5)
exc1_syn.w = ‘rand() * mV’ # Random weights
exc2_syn = Synapses(neurons, neurons, model=eqs_exc2_synapse, on_pre=‘I_exc_post += w’)
exc2_syn.connect(p=0.5)
exc2_syn.w = ‘rand() * mV’ # Random weights
inh_syn = Synapses(neurons, neurons, model=eqs_inh_synapse, on_pre=‘i_inh_post += w’)
inh_syn.connect(p=0.5)
inh_syn.w = ‘-rand() * mV’ # Random negative weights
M = StateMonitor(neurons, ‘I_exc’, record=True)
Hi @elnaz91,
Can you please post the actual code you are running? You can do this by just pasting your code inside triple backtick quotes (see example here)
But anyway, I think you’re looking for summed variables. Have a look at the documentation, but I think it would be something like
from brian2 import *
tau_membrane = 20 * ms
eqs_neurons = '''
dv/dt = (-v + I_exc + I_inh) / tau_membrane : volt
I_exc : volt
I_inh : volt
'''
tau_exc1 = 10 * ms
eqs_exc1_synapse = '''
dx/dt = -x / tau_exc1 : 1
I_exc_post = w * x : volt (summed)
w : volt
'''
eqs_exc2_synapse = '''
w : volt
'''
eqs_inh_synapse = '''
w : volt
'''
neurons = NeuronGroup(10, eqs_neurons, threshold='v > 0 * mV', reset='v = 0 * mV')
exc1_syn = Synapses(neurons, neurons, model=eqs_exc1_synapse, on_pre='x += 1')
exc1_syn.connect(p=0.5)
exc1_syn.w = 'rand() * mV' # Random weights
exc1_syn.x = 10
exc2_syn = Synapses(neurons, neurons, model=eqs_exc2_synapse, on_pre='I_exc_post += w')
exc2_syn.connect(p=0.5)
exc2_syn.w = 'rand() * mV' # Random weights
inh_syn = Synapses(neurons, neurons, model=eqs_inh_synapse, on_pre='I_inh_post += w')
inh_syn.connect(p=0.5)
inh_syn.w = '-rand() * mV' # Random negative weights
M = StateMonitor(neurons, 'I_exc', record=True)
spkmon = SpikeMonitor(neurons)
run(100 * ms)
figure()
plot(M.t, M.I_exc[0])
figure()
plot(spkmon.t, spkmon.i, '.')
show()
Hi @pabloabur , Thank you for the suggestion, but I can’t edit my reply (I believe because I am not a trusted user yet).
However, I have not solved my issue yet. My goal is to implement these membrane equations in Brian2; since I am new to both Brian2 and neuroscience, I am still struggling.
dv= (-v +WS(t) - W_{ei}S_i(t)+ W_{ee}h(t)*S_e(t)
In this model, the membrane potential V(t) is influenced by both excitatory and inhibitory inputs through distinct synaptic pathways:
Feedforward (FF) Excitation: Each spike in the feedforward connections contributes to an increase in V(t) by a weight W. This weight represents the strength of the synaptic connection from the feedforward excitatory population to the neuron.
Lateral Excitation: The excitatory-to-excitatory (E to E) lateral connections are modeled with a weight W_{ee} . The effect of these connections on the membrane potential is determined by convolving the spike train S_e(t) from the excitatory population with an exponential kernel h(t) , resulting in a modified spike train that influences V(t) through the W_{ee} weight.
Inhibitory Input: The membrane potential is also affected by inhibitory inputs where each spike from the inhibitory population S_i(t) decreases V(t) by a weight W_{ei}. This weight represents the synaptic strength from the inhibitory to the excitatory population.
To monitor the dynamics of these interactions, the model needs to track excitatory I_{exc} and inhibitory I_{inh} currents separately. However, there is a reported issue in the implementation where I_{exc} increases indefinitely. This unbounded increase suggests a potential issue with how excitatory currents are integrated over time or a lack of appropriate decay mechanisms in the model.
Hi @elnaz91 , thanks for the clarification. I think I understand your problem/model. Can you post your code though? As you said, it sounds like a problem with the integration of excitatory currents, but it’s hard to understand what’s going on without the code. Note that I_exc_post += w * x : volt and dx/dt = -x / tau_exc1 : 1 (event-driven) as you posted earlier lead to errors, so I changed them in my previous post. Maybe that helps.
cheers
@pabloabur Thank you so much. It worked without error. But how can I test it to make sure it works correctly? (I mean testing the logic)
by the way my s_model according to your previous response is as follows:
s_model = '''
dx_syn/dt = -x_syn / tau_slow : 1 (clock-driven)
w : volt
'''
on_pre = 'I_exc += w*x_syn'
You could connect a Poisson group to the network and use StateMonitors to record I_{exc}, I_{inh}, and v. Excitatory currents should look like a filtered version of inhibitory currents, and the membrane should evolve accordingly. Just try using low-frequency inputs or sparser connections so you don’t get lost when neurons start to spike. I think that’s a simple way to test it; does it make sense to you?
Yes, that makes perfect sense.
Thank you so much. It seems that works correctly.
and for the last question, can we interpret this synapse as weighted firing rates of a population?
I don’t think so. You could use a monitor to record the firing rate of a population though, as shown here
yes, that’s correct, but this function is for monitoring a population firing rate, not for using weighted firing rates for synaptic input.
Sorry for taking your time and I really appreciate your help.
Thanks.
I actually came up with the idea if adding two synapse and implement it as follows:
from brian2 import *
# Parameters
tau_slow = 10*ms # Time constant for slow exponential kernel
tau = 10*ms # Time constant for post-synaptic current decay
wei = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]) * 0.01 # Synaptic weight
wei = wei.flatten()
# Number of neurons
N = 3
rate = 100*Hz # Rate of the Poisson input
# Neuron group with excitatory neurons
exc = NeuronGroup(N, '''
dv/dt = (-v + I_exc) / tau : volt
dI_exc/dt = -I_exc / tau : volt
dx_syn/dt = -x_syn / tau_slow : volt
''',
threshold='v > 0.5*volt', reset='v = 0*volt', method='euler')
exc.v = 0*volt
exc.x_syn = 0*volt
# Add Poisson input to each neuron
poisson_input = PoissonInput(exc, 'v', N=1, rate=rate, weight=1*volt)
# Synapses with exponential kernel convolution
a = Synapses(exc, exc,
on_pre='''
x_syn += 1*volt
''',
method='euler')
a.connect('i == j')
# New synapses group for I_exc update
syn_exc = Synapses(exc, exc, model='''
w : 1
''',
on_pre='''
I_exc_post += w * x_syn
''',
method='euler')
syn_exc.connect()
syn_exc.w = wei
# Monitor for recording variables
M = StateMonitor(exc, ['v', 'I_exc', 'x_syn'], record=True)
spike_monitor = SpikeMonitor(exc)
# Run the simulation
run(20*ms)
What do you think about it?
Hi @elnaz91. I might have misunderstood, but I think this model is more complicated than necessary. If your excitatory current model consists of convolving the spike train with an exponential, then this is already what your x_syn variable does. This is somewhat similar to the model we discuss in Converting from integrated form to ODEs — Brian 2 2.6.0 documentation. If you think of something like this (note that w_{spike} is a short-hand here for: weight of the synapse that transmitted the spike)
Then this would be implemented in Brian by
And increasing I_{syn} by w_{spike} for each incoming spike.
So, adapting your code a bit (and removing the noise input and letting every neuron spike once):
from brian2 import *
# Parameters
tau_slow = 10*ms # Time constant for slow exponential kernel
tau = 10*ms # Time constant for post-synaptic current decay
wei = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]) * 0.01 # Synaptic weight
wei = wei.flatten()
# Number of neurons
N = 3
rate = 100*Hz # Rate of the Poisson input
# Neuron group with excitatory neurons
exc = NeuronGroup(N, '''
dv/dt = (-v + I_exc) / tau : volt
dI_exc/dt = -I_exc / tau_slow : volt
''',
threshold='v > 0.5*volt', reset='v = 0*volt', method='euler')
exc.v = 0*volt
# Add Poisson input to each neuron
# poisson_input = PoissonInput(exc, 'v', N=1, rate=rate, weight=1*volt)
# Synapses with exponential kernel convolution
syn_exc = Synapses(exc, exc,
"w : volt",
on_pre='''
I_exc += w
''',
method='euler')
syn_exc.connect()
syn_exc.w = wei*volt
# Monitor for recording variables
M = StateMonitor(exc, ['v', 'I_exc'], record=True)
spike_monitor = SpikeMonitor(exc)
# Run the simulation
run(20*ms)
exc.v[0] = 1*volt
run(20*ms)
exc.v[1] = 1*volt
run(20*ms)
exc.v[2] = 1*volt
run(40*ms)
will give me the following (the three colors represent the three neurons):
Is this what you wanted to achieve, or did I misunderstand your model?
Thank you, @mstimberg! That is totally correct, but the reason I split this into two synapses is that I have another synapse, which I did not include in this code. In that synapse, the I_{\text{post}} is increased by the amount of synaptic weight on the pre-synaptic spike. At the end of the experiment, I want to track all these excitatory incoming currents to each neuron in this population.
syn_fast = Synapses(input, exc, model='w:1', on_pre='I_exc += w*volt')
syn_fast.connect()
syn_fast.w = self.W4[0].flatten()
In general, the mathematics of my population differential equation is as follows:
\tau \frac{dv}{dt} = -v + \sum_{\text{preneuron}} w_{\text{slow}} h(t) \delta(t) + w_{\text{fast}} \delta(t)
The difference between fast and slow synapses is that in the fast synapse, the I_{exc} ccurrent is increased by the amount of w_{fast} when a pre-synaptic spike occurs. In the slow synapse, the I_exc current for each neuron at each time point increases by the amount of the weighted firing rates of all connected neurons. Since we have an all-to-all connection here, it would be the weighted firing rates of the entire population (clock-driven).
Hi again. Do you mean that the slow and fast synapses are triggered by the same spikes/connections? If yes, you could use a single Synapses object with two lines for on_pre. In either case, if I understand your model of w_{fast} correctly, it should be implemented as on_pre='v += w*volt'.
Hi, thank you for your prompt response. No, not the same connections. different connections, different weights.
About the fast one, if I use
on_pre = 'v+=w*volt
in that case, I could not track these currents to my neurons.
In fact, my goal is to study the excitatory incoming currents and, later, the inhibitory incoming currents to each neuron in the population. Right now, I don’t care if the I_exc is increased because of fast or slow currents. I just want to check how this current is changing over time.
Ahh, ok, now I get the issue. Sorry. Indeed, you cannot easily record this “current”, since it is an infinitely short event, not a current that has a value over the duration of the time step.
There is more than one possible workaround for this issue. Maybe the clearest I could think of would be something like this (ignoring the slow synapse):
eqs = """
dv/dt = -v/tau : 1 # does not refer directly to I_fast!
I_fast : 1 # total fast synaptic "current"
"""
neurons = NeuronGroup(...)
syn = Synapses(neurons, neurons, "w : 1", on_pre="I_fast += w")
# Update v and reset fast current for next timestep
neurons.run_regularly("""v += I_fast
I_fast = 0""", when="after_synapses", dt=defaultclock.dt)
If you want to record I_fast, you just have to be careful to record it at the right point during the simulation run (after it has been updated by the synapses, but before it has been reset). See Running a simulation — Brian 2 2.6.0 documentation for details. To make your life easier, you could also split the operation in two, so that the value I_fast is accessible by a StateMonitor at the start or end of a time step (the usual time for recordings):
# Reset previous current before applying synapses
neurons.run_regularly("I_fast = 0", when="before_synapses", dt=defaultclock.dt)
# Apply new currents after applying synapses
neurons.run_regularly("v += I_fast", when="after_synapses", dt=defaultclock.dt)
Let me know whether this makes sense or not 
Thank you!
Based on what we discussed, I wrote a code that is just an example of what I am doing. It includes two exc and inh populations and feedforward poisson input. All synapses are fast, and just EE synapses are slow. What do you think about it?
I actually prefer not to have this run_regurarly because I am implementing a hierarchical network, and this is just one layer of that.
from brian2 import *
# Parameters
tau_slow = 10*ms # Time constant for slow exponential kernel
tau = 10*ms # Time constant for post-synaptic current decay
# Number of neurons
N = 3
rate = 100*Hz # Rate of the Poisson input
eqs = '''
dv/dt = (-v + I_slow) / tau : volt
dI_slow/dt = -I_slow / tau : volt
I_exc = I_slow + I_fast : volt
I_fast : volt
I_inh : volt
'''
# Neuron group with excitatory neurons
exc = NeuronGroup(N, eqs, threshold='v > 0.5*volt', reset='v = 0*volt', method='euler')
exc.v = 0*volt
# Inhibitory population
inh = NeuronGroup(N, eqs, threshold='v > 0.5*volt', reset='v = 0*volt', method='euler')
inh.v = 0*volt
# Define the Poisson spike train for excitatory neurons
poisson_group = PoissonGroup(10, rates=rate)
# Connect PoissonGroup to excitatory neurons
poisson_syn = Synapses(poisson_group, exc, model='w:1', on_pre='I_fast += w*volt') # fast feedforward
poisson_syn.connect()
syn_ee = Synapses(exc, exc, "w : 1", on_pre="I_slow += w*volt") # excitatory ei synapse (slow)
syn_ee.connect()
syn_ee.w = np.random.rand(3,3).flatten()
syn_ie = Synapses(inh, exc, "w : 1", on_pre="I_inh -= w*volt") # Inhibitory ie synapse (fast)
syn_ie.connect()
syn_ie.w = np.random.rand(3,3).flatten()
syn_ei = Synapses(exc, inh, "w : 1", on_pre="I_fast += w*volt") # excitatory ei synapse (fast)
syn_ei.connect()
syn_ei.w = np.random.rand(3,3).flatten()
syn_ii = Synapses(inh, inh, "w : 1", on_pre="I_inh -= w*volt") # Inhibitory ii synapse (fast)
syn_ii.connect()
syn_ii.w = np.random.rand(3,3).flatten()
# Reset previous current before applying synapses
exc.run_regularly("I_fast = 0*volt; I_inh=0*volt", when="before_synapses", dt=defaultclock.dt)
# Apply new currents after applying synapses
exc.run_regularly("v += I_fast; v -= I_inh", when="after_synapses", dt=defaultclock.dt)
# Reset previous current before applying synapses
inh.run_regularly("I_fast = 0*volt; I_inh=0*volt", when="before_synapses", dt=defaultclock.dt)
# Apply new currents after applying synapses
inh.run_regularly("v += I_fast; v -= I_inh", when="after_synapses", dt=defaultclock.dt)
# Monitor for recording variables
M = StateMonitor(exc, ['v', 'I_exc'], record=True)
spike_monitor_exc = SpikeMonitor(exc)
spike_monitor_inh = SpikeMonitor(inh)
# Run the simulation
run(20*ms)
figure(figsize=(15, 6))
# Plot post-synaptic current
subplot(121)
plot(M.t/ms, M.I_exc[0], label='I_exc')
xlabel('Time (ms)')
ylabel('I_exc')
legend()
# Plot spike trains
subplot(122)
plot(spike_monitor_exc.t/ms, spike_monitor_exc.i, '.k', label='Excitatory neurons')
plot(spike_monitor_inh.t/ms, spike_monitor_inh.i, '.r', label='Inhibitory neurons')
xlabel('Time (ms)')
ylabel('Neuron index')
title('Spike trains')
legend()
show()
PS. code runs without error with no response.
Hi @elnaz91 . From a cursory glance, this looks reasonable to me. You are not getting any spikes, since you did not set the weights of poisson_syn (i.e., the weights are zero). But you should also be careful of how you connect things (or use PoissonInput instead of PoissonGroup): in your example, you fully connect the 10 Poisson neurons to the three excitatory neurons. This means, that each of the three neurons will get input from the same 10 neurons, i.e. all neurons will get the same external random input.
I don’t quite understand what you mean by that, but if you mean that you’d have to multiply these run_regularly operations for a network with multiple layers, then this is not the case. Instead of creating NeuronGroups/Synapses for every layer, you can create a single NeuronGroup for all layers, and use subgroups for the layers. In general, try to use the smallest number of objects. In your example, you could have a single NeuronGroup for both excitatory and inhibitory neurons:
neurons = NeuronGroup(2*N, eqs, threshold='v > 0.5*volt', reset='v = 0*volt', method='euler')
# Excitatory neurons
exc = neurons[:N]
# Inhibitory population
inh = neurons[N:]