Description of problem
I’m attempting to model a couple of neurons using Izhikevich 2003, with biexponential synapses for NMDA, AMPA, and GABA-A. I’ve included all the variables related to the synapses in the neuron model which I’m not sure is the correct approach. However, when I’ve tried to define the synaptic current through the synapse equations, I’ve not been able to get it to run and get different types of errors depending on how I define it.
I wanted to see if I have understood correctly how Brian defines the model and if this is an accurate way of representing the synapses – and if I can scale this up by simply defining more neurons within the population and later introduce intraconnections and connections to other neuron groups.
I also keep getting a warning log along the lines of ‘INFO:root:building ‘cython_magic’ extension’ and ‘INFO:root:clang’ and ‘INFO:root:creating’ which looks concerning but doesn’t seem to be an actual error. These warnings also sometimes don’t show at all, particularly when I re-run it after the first time. Should I be concerned?
Thank you in advance!
Minimal code to reproduce problem
def calculate_t_peak(tau_decay,tau_rise):
return ((tau_decay*tau_rise)/(tau_decay-tau_rise))*np.log(tau_decay/tau_rise)
def calculate_f(t_peak,tau_decay,tau_rise):
f = 1/(exp(-t_peak/tau_decay)-exp(-t_peak/tau_rise))
return f
defaultclock.dt = 0.01*ms
start_scope()
Ne = 1
Ni = 1
N = Ne + Ni
# AMPA- and NMDA-mediated neuron (excitatory)
tau_e_rise_ampa = 0.5*ms
tau_e_decay_ampa = 2.4*ms
t_peak = calculate_t_peak(tau_e_decay_ampa,tau_e_rise_ampa)
f_e_ampa = calculate_f(t_peak,tau_e_decay_ampa,tau_e_rise_ampa)
g_peak_e_ampa = 0.05 # synaptic weight
tau_e_rise_nmda = 8*ms
tau_e_decay_nmda = 200*ms
t_peak = calculate_t_peak(tau_e_decay_nmda,tau_e_rise_nmda)
f_e_nmda = calculate_f(t_peak,tau_e_decay_nmda,tau_e_rise_nmda)
g_peak_e_nmda = 0.075 # synaptic weight
# GABA-A-mediated neuron (inhibitory)
tau_i_rise = 1 * ms
tau_i_decay = 7 * ms
t_peak = calculate_t_peak(tau_i_decay,tau_i_rise)
g_peak_i = 0.33 # synaptic weight
f_i = calculate_f(t_peak,tau_i_decay,tau_i_rise)
#### define neuron model (Izhikevich, 2003) and biexponential synapse model ###
eqs = """dv/dt = (0.04*v**2 + 5*v + 140 - u + I_offset + I_syn + I_noise)/ms : 1
du/dt = (a*(b*v - u))/ms : 1
I_offset : 1
I_noise : 1
f : 1
a : 1
b : 1
c : 1
d : 1
ds_e1_ampa/dt = -s_e1_ampa/tau_e_decay_ampa : 1
ds_e2_ampa/dt = -s_e2_ampa/tau_e_rise_ampa : 1
s_e_ampa = s_e1_ampa - s_e2_ampa : 1
ds_e1_nmda/dt = -s_e1_nmda/tau_e_decay_nmda : 1
ds_e2_nmda/dt = -s_e2_nmda/tau_e_rise_nmda : 1
s_e_nmda = s_e1_nmda - s_e2_nmda : 1
ds_i1/dt = -s_i1/tau_i_decay : 1
ds_i2/dt = -s_i2/tau_i_rise : 1
s_i = s_i1 - s_i2 : 1
g_nmda = 1 / (1 + exp(-v/16.13) / 3.57) : 1
I_syn_e_ampa = g_peak_e_ampa * f_e_ampa * s_e_ampa : 1
I_syn_e_nmda = g_nmda * g_peak_e_nmda*f_e_nmda*s_e_nmda : 1
I_syn_i = -g_peak_i*f_i*s_i : 1
I_syn = I_syn_e_ampa + I_syn_e_nmda + I_syn_i : 1
"""
reset_eq = 'v=c; u+= d;'
threshold = 'v>30'
P = PoissonGroup(Ne,rates=10*Hz)
G = NeuronGroup(N,eqs,method='euler',threshold=threshold,reset=reset_eq)
G_exc = G[:Ne]
G_inh = G[Ne:]
# Initialise variables
# Excitatory
G_exc.v = -70
G_exc.a = 0.02
G_exc.b = 0.2
G_exc.c = -65
G_exc.d = 2
G_exc.u = 'b*v'
G_exc.I_offset = 3.47
# Inhibitory
G_inh.v = -70
G_inh.a = 0.1
G_inh.b = 0.2
G_inh.c = -65
G_inh.d = 2
G_inh.u = 'b*v'
G_inh.I_offset = 3.85
S_e = Synapses(P,G_exc,on_pre='s_e1_ampa+=g_peak_e_ampa; s_e1_nmda +=g_peak_e_nmda')
S_i = Synapses(P,G_inh,on_pre='s_i1+=g_peak_i')
S_e.connect(j='i')
S_i.connect(j='i')
G_exc.run_regularly("I_noise = 5*randn()", dt=0.01 * ms)
G_inh.run_regularly("I_noise = 2*randn()", dt=0.01 * ms)
run(1000*ms)