Description of problem
I am trying to connect two neurons using alpha synapse for Izhikevich model (which I want to keep unitless). The code seems to be working but I am not getting the right plots. I don’t expect my code to be all right, there must be some issues. Any insights are much appreciated
Minimal code to reproduce problem
from brian2 import *
start_scope()
defaultclock.dt = 0.01*ms
num_neurons = 2
duration = 2*second
#Izikevich Model Parameters
a = 0.02
b = 0.2
c = -65
d = 8
I = 15
#Parameters
El = -60*mV
gl = 0.7*msiemens/cm**2
#Alpha Synapse
tau = 0.3*ms
g_synpk = 0.5
g_synmaxval = (g_synpk)/(tau*exp(-1)/ms)
eqs_il = '''
il = gl * ((El/volt)-v) :amp/meter**2/volt
'''
eqs = '''
dv/dt = (0.04*v**2 + 5*v + 140 -u + (I/amp) + GAMPA*(10-v))/ms : 1
du/dt = a*(b*v-u)/ms : 1
'''
reset = '''
v = c
u = u + d
'''
eqs_syn = '''
I : amp
GAMPA : 1
'''
eqs = eqs+eqs_il+eqs_syn
G = NeuronGroup(num_neurons,eqs,threshold='v>=30', reset=reset, method='euler')
G.v = c
G.u = b*c
Sr = Synapses(G, G, clock=G.clock, model='''
g_synmax:1
dg/dt = -g/tau + z : 1
dz/dt = -z/tau : Hz
GAMPA_post = g : 1 (summed)
''',
on_pre = '''
z +=g_synmax*Hz
'''
)
Sr.connect(i=[0],j=[1])
Sr.g_synmax = g_synmaxval
Sr.delay = 1*ms
M = StateMonitor(G,('v','GAMPA'), record=True)
G.I[0] = 0*nA
G.I[1] = 0*nA
run(5.0*ms,report='text')
G.I[0] = 8*nA
G.I[1] = 0*nA
run(5.0*ms,report='text')
G.I[0] = 0*nA
G.I[1] = 0*nA
figure
subplot(5,1,1)
plot(M.t/ms, M.v[0])
subplot(5,1,2)
plot(M.t/ms, M.v[0]*g_synpk)
subplot(5,1,3)
plot(M.t/ms, M.v[1]*g_synpk,'r')
subplot(5,1,4)
plot(M.t/ms, M.v[1],'r')
subplot(5,1,5)
plot(M.t/ms, M.GAMPA[1])
What you have aready tried
These are the plots I am getting.
Hi @shehnaz.
I don’t think there’s anything fundamentally wrong with the equations. One little note, though: if all synapses have the same values for tau, it would be more efficient to move the equations for the alpha synapse into the neuron instead of having them in the Synapses equations. In a network this can make a big difference, since the equations then have to be solved only once per every neuron, not once for every synapse. In your current example (with a single synapse), this does not matter, though.
You are not seeing any effect of your input current, since it is much too weak: you include I/amp in the equations to make the current dimensionless, but this means that your 8nA current is added as 8\cdot 10^{-9} to the equations, i.e. as a tiny value. With I/nA you should get something visible. The same is true for your g_synmaxval value, it needs to be much bigger to have a visible effect on the target cell (e.g. try multiplying it by 100).
Oh, and note that this code:
G.I[0] = 0*nA
G.I[1] = 0*nA
run(5.0*ms,report='text')
G.I[0] = 8*nA
G.I[1] = 0*nA
run(5.0*ms,report='text')
G.I[0] = 0*nA
G.I[1] = 0*nA
means that you run the simulation without input for 5ms, and then with input for 5ms. You probably intended to run it for 5ms again without input current, so this needs another run(...) in the end.
Two additionalminor comments:
- For better readability, please surround code examples by triple backticks like this:
```
# Your python code here
```
This way the code gets displayed nicely and the special characters don’t mess with the markdown syntax (I edited your post to do this).
- Your equations contain the definition of a leak current
il, but this current is not used anywhere in your membrane potential equations (it would be redundant with Izhikevich’s model any way) – you can simply remove it.
from brian2 import *
start_scope()
defaultclock.dt = 0.01*ms
num_neurons = 3
duration = 2*second
#Izikevich Model Parameters
a = 0.02
b = 0.25
c = -65
d = 2
I = 15
#Alpha Synapse
tau = 0.3*ms
g_synpk = 0.08
g_synmaxval = (g_synpk)/(tau*exp(-1)/ms)
eqs = '''
dv/dt = (0.04*v**2 + 5*v + 140 -u + (I/nA) + GAMPA*(10-v))/ms : 1
du/dt = a*(b*v-u)/ms : 1
I : amp
GAMPA : 1
'''
reset = '''
v = c
u = u + d
'''
G = NeuronGroup(num_neurons,eqs,threshold='v>=30', reset=reset, method='euler')
G.v = c
G.u = b*c
Sr = Synapses(G, G, clock=G.clock, model='''
g_synmax:1
dg/dt = -g/tau + (z/ms) : 1
dz/dt = -z/tau : 1
GAMPA_post = g : 1 (summed)
''',
on_pre = '''
z +=g_synmax
'''
)
Sr.connect(i=[0],j=[1])
Sr.g_synmax = g_synmaxval
Sr.delay = 1*ms
M = StateMonitor(G,('v','GAMPA'), record=True)
G.I[0] = 0*nA
G.I[1] = 0*nA
run(5.0*ms,report='text')
G.I[0] = 80*nA
G.I[1] = 0*nA
run(20.0*ms,report='text')
G.I[0] = 0*nA
G.I[1] = 0*nA
run(5.0*ms, report='text')
figure
subplot(6,1,1)
plot(M.t/ms, M.v[0])
subplot(6,1,2)
plot(M.t/ms, M.v[0]*g_synpk)
subplot(6,1,3)
plot(M.t/ms, M.v[1]*g_synpk,'r')
subplot(6,1,4)
plot(M.t/ms, M.v[1],'r')
subplot(6,1,5)
plot(M.t/ms, M.GAMPA[0])
subplot(6,1,6)
plot(M.t/ms, M.GAMPA[1])
I made the changes in the current and I modified my z to be dimensionless, now I am getting voltage spikes and current plot.
I want to connect the neuron 0 to neuron 1 using an inhibitory connection, and I am not sure how to do that.
This is the output I am getting after making the changes.
Thank you very much for your help and time. I really appreciate the way you explained the problems and small nuances.
I want to learn coding in brian2 but I am new to both neuroscience as well as coding. Any good resources that you could suggest would be a great help.
