Description of problem
I have implemented pinsky and rinzel neuron model in Brian and I am creating a microcircuit model wherein I connect a Wang and Buzsaki interneuron model to the pyramidal cell. I have two issues that I am not able to figure out.
- In the individual implementation of Pinsky and Rinzel cell, I dont see spiking and bursting behavior in the same simulation by manipulating the dendritic current. Bursts are only seen when condutance (gc) is increased to 4 which is higher than the value at which we orignally should see bursts (ie at 3mS). And I am not able to figure out why that is happening.
- In the implementation of the circuit, my synaptic gating variable steps out of its bound (it should be between 0 and 1, it goes upto to 1.6). I am guessing that there is some issue in the way the synapse equations are written but I cannot pinpoint the origin of the issue.
Here is a ref for better understanding of my circuit model, and for Pyramidal cell model:
- Modulation of spike and burst rate in a minimal neuronal circuit with
feed-forward inhibition (Fleur Zeldenrust, Wytse J. Wadman, 2012)
Minimal code to reproduce problem
Pinsky and Rinzel implementation:
#this is the edited version of the code Fleur forwarded on slack
#also the final one that seems to work - and no the orignal pinksy and rinzel paramter values - its the Fleur's version and the booth and bose version
start_scope()
pyr_parameters = {
"gLs" : 0.1*mS,
"gLd" : 0.1*mS,
"gNa" : 30.*mS,
"gKdr" : 15.*mS,
"gCa" : 10.*mS,
"gKahp" : 0.8*mS,
"gKC" : 15.*mS,
"VNa" : 60.0*mV,
"VCa" : 80.0*mV,
# "VCa" : 100.0*mV,
"VK" : -75.0*mV,
"VL" : -60.0*mV,
"Vsyn" : 0.0*mV,
"Cm" : 3.0*uF,
"gc" : 2.1*mS,
# "gc" : 2.1*mS, #for fig. 7: 3.0*mS, normal: 2.1*mS, for fast spiking 1.425*mS
"pp" : 0.5,
"Id" : -0.3*uA,
"Is" : 0.0*uA,
"alphac": 2.0/ms,
"betac" : 0.1/ms,
"betaqd": 0.001/ms,
}
pyr_initialvalues = {
'Vs' : -62.89223689*mV,
# 'Vs' : 35.0*mV,
# 'Vd' : 35.0*mV,
'Vd' : -62.98248752*mV,
# 'Vs' : 0*mV,
# 'Vd' : 0*mV,
'Cad' :0.21664282,
'hs' : 0.99806345,
'ns' :0.00068604,
'sd' : 0.01086703,
'cd' :0.00809387,
'qd' : 0.0811213,
}
eqs_pyr = '''
#membrane potential equations
#for soma
dVs/dt=(-ILs-INa-IK_DR+(gc/pp)*(Vd-Vs)+(Is/pp))/Cm: volt
#for dendrite
dVd/dt=(-ILd-ICad-IK_ahp-IK_C+((gc/(1.0-pp))*(Vs-Vd))+Id_FN/(1.0-pp))/Cm: volt #-Isynd/(1.0-pp)
ILs = gLs*(Vs-VL) : amp
INa = gNa*Minfs*Minfs*hs*(Vs-VNa) :amp
IK_DR = gKdr*ns*(Vs-VK) :amp
ILd = gLd*(Vd-VL) :amp
ICad = gCa*sd*sd*(Vd-VCa) :amp
IK_ahp = gKahp*qd*(Vd-VK) :amp
IK_C = gKC*cd*chid*(Vd-VK) :amp
#calcium concentration
dCad/dt= (-0.13*ICad/uA-0.075*Cad)/ms: 1
# soma dynamics
dhs/dt= alphahs-(alphahs+betahs)*hs: 1
dns/dt= alphans-(alphans+betans)*ns: 1
# dendrite dynamics
dcd/dt= alphacd-(alphacd+betacd)*cd: 1
dqd/dt= alphaqd-(alphaqd+betaqd)*qd: 1
# calcium dynamics
dsd/dt= alphasd-(alphasd+betasd)*sd: 1
# soma state
# sodium current
alphams= 0.32*(-46.9-Vs/mV)/(exp((-46.9-Vs/mV)/4.0)-1.0)/ms: Hz
betams= 0.28*(Vs/mV+19.9)/(exp((Vs/mV+19.9)/5.0)-1.0)/ms: Hz #orignal one
# betams= 0.28*(Vs/mV+70)/(exp((Vs/mV+70)/5.0)-1.0)/ms: Hz
Minfs= alphams/(alphams+betams): 1
alphahs= 0.128*exp((-43.0-Vs/mV)/18.0)/ms: Hz
betahs= 4.0/(1.0+exp((-20.0-Vs/mV)/5.0))/ms: Hz
# Kdr potassium current
alphans= 0.016*(-24.9-Vs/mV)/(exp((-24.9-Vs/mV)/5.0)-1.0)/ms: Hz
betans= 0.25*exp(-1.0-0.025*Vs/mV)/ms: Hz
# dendrite state
# calcium induced potassium current
alphacd=((int(Vd/mV<=-10)*(exp(((Vd/mV+50.0)/11) - (-53.3-Vd/mV) / 27))/18.975) + (int(Vd/mV>-10)*2.0*exp((-53.5-Vd/mV)/27.0)))/ms : Hz
betacd=(int(Vd/mV<=-10)*(2.0*(exp((-53.5-Vd/mV)/27.0))-alphacd*ms))/ms : Hz
# AHP: calcium induced potassium current
alphaqd= clip(0.00002*Cad,0.00002*Cad, 0.01)/ms: Hz
betaqd = 0.001/ms: Hz
chid= clip(Cad/250.,Cad/250. , 1. ) : 1
# calcium current
alphasd= (1.6/(1.0+exp(-0.072*(Vd/mV-5.0))))/ms: Hz
betasd= (0.02*(Vd/mV+8.9)/(exp((Vd/mV+8.9)/5.0)-1.0))/ms: Hz
# betasd= (0.02*(Vd/mV+13)/(exp((Vd/mV+13)/5.0)-1.0))/ms: Hz
# Is_FN = Iapp(t) :amp
Is : amp
Id : amp
Id_FN = Iapp(t) : amp
'''
#neuron groups
pyr_model = NeuronGroup(1, model=eqs_pyr, threshold='Minfs > 0.5', refractory=2 * ms,
reset=None, dt=dt,method='rk4',name='pyr',
namespace=pyr_parameters)
#initial values
pyr_model.set_states(pyr_initialvalues)
pyr_model.Vs = pyr_initialvalues["Vs"]
pyr_model.Vd = pyr_initialvalues["Vd"]
mon_pyr = StateMonitor(pyr_model, ['Vs','Vd','Cad','ILs',"INa",'ICad','IK_ahp','IK_C','ILd'], record=True)
I = StateMonitor(pyr_model, 'Id', record=True)
run(simtime)
fig,ax = plt.subplots(figsize=[12,4])
ax.plot(mon_pyr.t/ms, mon_pyr.Vs[0]/mV, label='Vs')
ax.plot(mon_pyr.t/ms, mon_pyr.Vd[0]/mV, label='Vd')
ax.set_xlabel('Time (ms)',fontsize=18)
ax.set_ylabel('Memb.Potential of Pyr (mV)',fontsize=18)
# ax1 = ax.twinx()
# ax1.plot(mon_pyr.t/ms,mon_pyr.Cad[0],c='g',label='Cad')
ax.legend(loc=[0.2,1])
ax.legend(loc=[0.1,1])
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.title('Pyramidal Cell Spike (frozen noise)',fontsize=18)
show()
Implementation of Synapses:
eqs_synapse = '''
Isyn_syn = gsyn * ssyn * (V_postsynapse - vsyn) : amp
dssyn/dt = alpha_syn * (1 - ssyn) / (1 + exp(-(V_presynapse - theta_syn) / sigma_syn)) - beta_syn * ssyn : 1 (event-driven)
alpha_syn = 1 / tau_rise : Hz
beta_syn = 1 / tau_decay : Hz
gsyn : siemens
vsyn : volt
tau_rise : second
tau_decay : second
theta_syn : volt
sigma_syn : volt
V_postsynapse : volt
V_presynapse : volt
'''
synapse_ff = Synapses(int_model, pyr_model, model=eqs_synapse, on_pre="Id += Isyn_syn" ,method='euler',dt=dt)
synapse_ff.connect(i=0, j=0)
synapse_ff.gsyn = 8.0*mS
synapse_ff.vsyn = -62.0*mV
synapse_ff.tau_rise = 5.0*ms
synapse_ff.tau_decay = 20.0*ms
synapse_ff.theta_syn = 0.0*mV
synapse_ff.sigma_syn = 1.0*mV
synapse_ff.V_postsynapse = pyr_model.Vs
synapse_ff.V_presynapse = int_model.v_int
'''
(all equations used are mentioned in the reference mentioned above)
What you have aready tried
Expected output (if relevant)
Actual output (if relevant)
(for the pyr plots, Id = -0.3 , a value at which the cell should burst, but instead it spikes (gc = 2.1))
plot of the complete circuit in feedforward mode:
