Problem creating Pinsky and Rinzel Cell

Support
1

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,
    'Cad' : 0.,
    '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
dVd/dt=(-gLd*(Vd-VL)-ICad-gKahp*qd*(Vd-VK)-gKC*cd*chid*(Vd-VK)+(gc*(Vs-Vd)+Id+Isyn)/(1.0-pp))/Cm: volt

#calcium concentration
dCad/dt=  (-0.13*ICad/amp-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
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
alphaqd=  clip(0.00002*Cad, 0, 0.01)/ms: Hz
chid=     clip(Cad/250.0, 0, 1.0): 1

# calcium current
ICad=     gCa*sd*sd*(Vd-VCa): amp
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()

What you have aready tried

Expected output (if relevant)

image
image777×785 24 KB

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

Actual output (if relevant)

image
image1008×371 13.1 KB

Full traceback of error (if relevant)

2

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

which looks reasonable to me.

I am also not sure about this line:

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 :slight_smile: