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:

3

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,
    'Cad' :0.0,
    '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/(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/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 #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
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

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"]

mon_pyr = StateMonitor(pyr_model, ['Vs','Vd','Cad','ILs',"INa",'ICad','IK_ahp','IK_C','ILd'], record=True)

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.ICad[0]/mV, label='ICad')
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.plot(mon_pyr.t/ms,mon_pyr.Cad[0],c='g',label='Cad')
ax1.legend(loc=[0.2,1])
ax.legend(loc=[0.1,1])
show()

image
image1021×371 25.9 KB

image
image1075×417 45.5 KB

4

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:

Figure_1
Figure_11200×400 31.9 KB

Figure_2
Figure_21200×400 48.6 KB