# Description of problem

Hello, I am trying to create a Pinsky and Rinzel pyramidal cell model (the entire code for the same is pasted below). I am not able to figure out why the model does not show any spiking or bursting behavior even after providing substantial external input. I have rechecked my equations, parameter values and units however I am still unable to figure out the issue, maybe someone here could help me to understand why the cell does not show any activity.
Any help would be appreciated, thank you!

Reference papers for equations:

1. https://doi.org/10.1016/j.neunet.2012.12.008
2. Intrinsic and Network Rhythmogenesis in a Reduced Traub Model
for CA3 Neurons (Pinsky and Rinzel, 1994).

# Minimal code to reproduce problem

``````from brian2 import *

pyr_parameters = {
"Cm" : 3.*pF,
"gc" : 2.1*nS,
"pp" : 0.5,
"Is" : 0.*namp,
"Id" : 0.*namp,
"gLs": 0.1*nS,
"gLd": 0.1*nS,
"gNa": 30.*nS,
"gKdr":15.*nS,
"gCa": 10.*nS,
"gKahp": 0.8*nS,
"gKC": 15.*nS,
"VNa": 60.*mV,
"VCa": 80.*mV,
"VK" : -75.*mV,
"VL" : -60.*mV,
"Vsyn" : 0.*mV,
"alphac": 2./ms,
"betac" : 0.1/ms,
"betaqd": 0.001/ms,
"dt" : defaultclock.dt,
"simtime": 1*second,
"Ibaseline": 0.*namp,
"Iscale": 0.*namp,
"Isyn": 30.0*uA,
}
pyr_initialvalues = {
'Vs' : -70.*mV,
'Vd' : -70.*mV,
'hs' : 1.,
'ns' : 0.5,
'sd' : 0.,
'cd' : 0.5,
'qd' : 0.
}

Cm  = pyr_parameters['Cm']
gc  = pyr_parameters['gc']
pp  = pyr_parameters['pp']
Is  = pyr_parameters['Is']
Id  = pyr_parameters['Id']
gLs = pyr_parameters['gLs']
gLd  = pyr_parameters['gLd']
gNa  = pyr_parameters['gNa']
gKdr  = pyr_parameters['gKdr']
gCa  = pyr_parameters['gCa']
gKahp  = pyr_parameters['gKahp']
gKC  = pyr_parameters['gKC']
VNa  = pyr_parameters['VNa']
VCa  = pyr_parameters['VCa']
VK  = pyr_parameters['VK']
VL  = pyr_parameters['VL']
Vsyn  = pyr_parameters['Vsyn']
alphac  = pyr_parameters['alphac']
betac  = pyr_parameters['betac']
betaqd  = pyr_parameters['betaqd']

eqs_pyr = '''
#membrane potential equations
#for soma
dVs/dt=(-gLs*(Vs-VL)-gNa*Minfs*Minfs*hs*(Vs-VNa)-gKdr*ns*(Vs-VK)+(gc/pp)*(Vd-Vs)+Is/pp)/Cm: volt

#for dendrite

#calcium concentration

# 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
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)-(Vd/mV+53.5)/27)/18.975+int(1.0*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

# 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

Isyn : amp
'''
#time settings
defaultclock.dt = pyr_parameters["dt"]
duration = pyr_parameters["simtime"]

#neuron groups
pyr_model = NeuronGroup(1, model=eqs_pyr, threshold='Minfs > 0.5', refractory=2 * ms, reset=None, dt=pyr_parameters["dt"],method='euler',name='pyr')
pyr_model.Vs = pyr_initialvalues["Vs"]
pyr_model.Vd = pyr_initialvalues["Vd"]

mon_pyr = StateMonitor(pyr_model, ['Vs','Vd','Isyn'], record=True)

run(duration)

#plot the results
figure(figsize=(12,4))
plot(mon_pyr.t/ms, mon_pyr.Vs[0]/mV, label='Vs')
plot(mon_pyr.t/ms, mon_pyr.Vd[0]/mV, label='Vd')
xlabel('Time (ms)')
ylabel('v (mV)')
legend()
show()
``````

# Expected output (if relevant)

This image is from the Pinsky and Rinzel which I am trying to replicate.

# Full traceback of error (if relevant)

Hi @bindok. Complex models like this are tricky to debug! I’d recommend verifying them component by component and see whether all components are as you expected. One useful technique is described here: How to plot functions — Brian 2 2.6.0 documentation

In your code, it shows that some of the rate functions are incorrect, e.g. if you use the following instead of running your model:

``````pyr_model = NeuronGroup(1000, model=eqs_pyr, threshold='Minfs > 0.5', refractory=2 * ms, reset=None, dt=pyr_parameters["dt"],method='euler',name='pyr')
pyr_model.Vs = pyr_initialvalues["Vs"]
pyr_model.Vd = pyr_initialvalues["Vd"]

pyr_model.Vs = np.linspace(-80, 60, len(pyr_model))*mV
plt.plot(pyr_model.Vs/mV, pyr_model.alphams, label='alphams')
plt.plot(pyr_model.Vs/mV, pyr_model.betams, label='betams')
plt.legend()
plt.show()
``````

you will get

This is clearly not correct, since rate functions should always be positive.

Here,

``````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
``````

should be

``````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
``````

And similarly for the other rate functions of this form.

This gives you

which looks reasonable to me.

``````dCad/dt=  (-0.13*ICad/amp-0.075*Cad)/ms: 1
``````

I think you might need to use `pA` instead of `amp`?

A minor point (that doesn’t affect the results): I understand that you want to keep the parameters together in a dictionary. In that case, you don’t actually need to unpack them into variables of the same name as you do here:

``````Cm  = pyr_parameters['Cm']
gc  = pyr_parameters['gc']
pp  = pyr_parameters['pp']
...
``````

Instead, you can directly hand over the dictionary to the `NeuronGroup` as the `namespace` argument:

``````pyr_model = NeuronGroup(1, model=eqs_pyr, threshold='Minfs > 0.5', refractory=2 * ms,
reset=None, dt=pyr_parameters["dt"],method='euler',name='pyr',
namespace=pyr_parameters)
``````

Similarly, you can initialize many variables at once:

``````pyr_model.set_states(pyr_initialvalues)
``````

Hope that helps you get unstuck! Final remark: if you get to reproduce some representative results from the Pinsky and Rinzel paper, the code would be a nice addition to the Brian examples

I apologise for the late response. However, I was able to make the changes in order to get the spiking activity. But the issue currently is the exponential increase of the currents at the end of the simulation. I believe this is a common issue with HH type neuron models so I was wondering how are they usually resolved?
Here is my new updated code

``````from brian2 import *
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" : 55.0*mV,
"VCa" : 80.0*mV,
# "VCa" : 100.0*mV,
"VK" : -75.0*mV,
"VL" : -60.0*mV,
"Vsyn" : 0.0*mV,
"Cm" : 3.*uF,
"gc" : 1.425*mS,
# "gc" : 2.1*mS,
"pp" : 0.5,
"Id" : 0.0*uA,
"Is" : 0.5*uA,
# "Id" : 0.0*uA,
# "Is" : -0.3*uA, / -0.5*uA
"alphac": 2.0/ms,
"betac" : 0.1/ms,
"betaqd": 0.001/ms,
"dt" : 0.01*ms,
"simtime": 300*ms,
# "Ibaseline": 0.*uA,
# "Iscale": 0.*uA,
"Isyns": 2.0*uA,
"Isynd": 0.0*uA,
}

pyr_initialvalues = {
'Vs' : -62.89223689*mV,
'Vd' : -62.98248752*mV,
# 'Vs' : -4.6*mV,
# 'Vd' : -4.5*mV,
'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

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
IK_ahp = gKahp*qd*(Vd-VK) :amp
IK_C   = gKC*cd*chid*(Vd-VK) :amp

#calcium concentration

# 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 #fuck this

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<=-50)*((exp((Vd/mV-10.0)/11))-(exp(Vd/mV-6.5)/27))/18.97+int(1.0*Vd/mV>-50)*2.0*exp((6.5-Vd/mV)/27.0))/ms  : Hz
betacd=(int(Vd/mV<=-50)*(2.0*(exp((6.5-Vd/mV)/27.0))-alphacd*ms))/ms : Hz

# AHP: calcium induced potassium current
betaqd = 0.001/ms: Hz

# 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

Isyns : amp
Isynd : amp
Iext : amp
'''

#time settings
defaultclock.dt = pyr_parameters["dt"]
duration = pyr_parameters["simtime"]

#neuron groups
pyr_model = NeuronGroup(1, model=eqs_pyr, threshold='Minfs > 0.5', refractory=2 * ms,
reset=None, dt=pyr_parameters["dt"],method='euler',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"]

run(duration)

#plotting ICad, ILd, IK_ahp, IK_C, ILd
figure(figsize=(12, 4))
plot(mon_pyr.t/ms, mon_pyr.ILd[0]/mV, label='ILd')
plot(mon_pyr.t/ms, mon_pyr.IK_ahp[0]/mV, label='IK_ahp')
plot(mon_pyr.t/ms, mon_pyr.IK_C[0]/mV, label='IK_C')
plot(mon_pyr.t/ms, mon_pyr.ILs[0]/mV, label='ILs')
plot(mon_pyr.t/ms, mon_pyr.ILd[0]/mV, label='ILd')
legend()
xlabel('Time (ms)')
ylabel('Current (mV)')
show()

#plot the results
fig,ax = plt.subplots(figsize=[12,4])
ax.plot(mon_pyr.t/ms, -60+mon_pyr.Vs[0]/mV, label='Vs')
ax.plot(mon_pyr.t/ms, -60+mon_pyr.Vd[0]/mV, label='Vd')

ax.set_xlabel('Time (ms)')
ax.set_ylabel('v (mV)')
ax1 = ax.twinx()
ax1.legend(loc=[0.2,1])
ax.legend(loc=[0.1,1])
show()
``````

Hi @bindok. I don’t think this is a general issue, but rather a problem with the `ICad`/`Cad` definition. In your equation, you have

``````dCad/dt=  (-0.13*ICad/amp-0.075*Cad)/ms: 1
``````

but I don’t think this is correct, since the original equations did not consider the `ICad` values to be in Ampère. With the conversions you did elsewhere, shouldn’t it rather be

``````dCad/dt=  (-0.13*ICad/uA-0.075*Cad)/ms: 1
``````

I’m not quite sure, but this gives the following which at least does not look completely unreasonable to me: