4
$\begingroup$

Not sure if this is the right place to ask this but I've been attempting to implement the Pinsky-Rinzel model in the Brian2 simulator. I'm basing my implementation off the original article here: http://www.caam.rice.edu/~cox/neuro/pinskyrinzel.pdf and I've borrowed some of the formalism from the XPP version of the model here: http://senselab.med.yale.edu/ModelDB/showmodel.cshtml?model=35358&file=/b04feb12/booth_bose.ode (Though they are shifted 60mV from where I'm working)

I've got it to replicate the somatic firing (as you can see if you set gc to zero), but I can't figure out how to get the complex dendritc and somatic interplay I'm supposed to see as seen in figure 2 and figure 3 of the model and my q and Ca values don't replicate the model at all. I think there is something wrong in how I'm handling the calcium currents, or the units I'm using there.

Does anyone have any ideas of where I'm going wrong?

My code:

# -*- coding: utf-8 -*-
"""
Created on Sat Oct 24 11:21:30 2015

@author: Conor
"""

from __future__ import division
import numpy as np
from brian2 import *
defaultclock.dt = 0.05*ms

gbarl=0.1*msiemens
gbarna=30.0*msiemens
gbarkDr=15.0*msiemens
gbarCa=10.0*msiemens
gbarkc=15.0*msiemens

gkahp=0.8*msiemens
gNMDA=0*msiemens
gampa=0*msiemens

Vna=120.0*mV
Vca=140.0*mV
Vk=-15.0*mV
Vl=0*mV
Vsyn=60*mV
Id=0*uamp
Is=.15*uamp
gc=2.1*msiemens
p=0.5
Cm=3*uF
eqs='''
alpham=0.32*(13.1-Vs/mV)/(exp((13.1-Vs/mV)/4)-1)/ms : Hz
betam=0.28*(Vs/mV-40.1)/(exp((Vs/mV-40.1)/5)-1)/ms : Hz
alphan=0.016*(35.1-Vs/mV)/(exp((35.1-Vs/mV)/5)-1)/ms : Hz
betan=0.25*exp(0.5-0.025*Vs/mV)/ms : Hz
alphah=0.128*exp((17-Vs/mV)/18)/ms : Hz
betah=4/(1+exp((40-Vs/mV)/5))/ms : Hz
alphas=1.6/(1+exp(-0.072*(Vd/mV-65)))/ms : Hz
betas=.02*(Vd/mV-51.1)/(exp((Vd/mV-51.1)/5)-1)/ms : Hz
alphac=(exp((Vd/mV-10)/11)-exp((Vd/mV-6.5)/27))/(18.975)* (Vd<=50*mV)/ms+2*exp((6.5-Vd/mV)/27)*(Vd>50*mV)/ms : Hz
betac=((2*exp((6.5-Vd/mV)/27))/ms-alphac)*(Vd<(50*mV)) +0/ms : Hz
alphaq=clip(.00002*Ca,0,0.01)/ms :Hz
betaq=0.001/ms : Hz
Xca=clip(Ca/250,0,1) : 1
dh/dt= (alphah*(1-h)-betah*h) : 1
dn/dt = (alphan*(1-n)-betan*n) : 1
ds/dt = alphas*(1-s)-betas*s : 1
dc/dt = alphac*(1-c)-betac*c : 1
dq/dt= (alphaq*(1-q)-betaq*q) : 1
dCa/dt=-0.13*Ica/uamp/ms-0.075*Ca/ms : 1
Minfs=alpham/(alpham+betam) : 1
ILeakVs=gbarl*(Vs-Vl) : amp
ILeakVd=gbarl*(Vd-Vl) : amp
INa=gbarna*(Minfs**2)*h*(Vs-Vna) : amp
Ik_dr=gbarkDr*n*(Vs-Vk) : amp
Ica=gbarCa*s*s*(Vd-Vca) : amp
IK_c=gbarkc*Xca*c*(Vd-Vk) : amp
IK_AHP=gkahp*q*(Vd-Vk) : amp
dVs/dt=(-ILeakVs-INa-Ik_dr+gc/p*(Vd-Vs)+Is/p)/Cm : volt
dVd/dt=(-ILeakVd-Ica-IK_AHP-IK_c+gc/(1-p)*(Vs-Vd)+INMDA/(1-p))/Cm : volt
INMDA=gNMDA*Si*(1+0.28*exp(-0.062*(Vs/mV-60)))**(-1)*(Vd-Vsyn) : amp
Hxs=0<=((Vs/mV-10)*1) : 1
Hxw=0<=((Vs/mV-20)*1) : 1
dSi/dt=(Hxs-Si/150)/ms : 1
'''
neurons = NeuronGroup(1, eqs)
neurons.Vs=-4.6*mV
neurons.Vd=-4.7*mV
neurons.h=.999
neurons.n=.0001
neurons.s=.009
neurons.c=.007
neurons.q=.01
neurons.Ca=0
neurons.Si=.3                            
M = StateMonitor(neurons, True ,record=True,dt=.1*ms)

run(500*ms, report='text') 
plt.cla()
plt.clf()
plot(M.t/ms, M[0].Vs/mV)
plot(M.t/ms, M[0].Vd/mV)
show()
$\endgroup$
3
$\begingroup$

If anyone ever finds this and has the same problem: The original paper has a flaw in the paper:

alphac=(exp((Vd/mV-10)/11)-exp((Vd/mV-6.5)/27))/(18.975)* (Vd<=50*mV)/ms+2*exp((6.5-Vd/mV)/27)*(Vd>50*mV)/ms : Hz

An errata published later changes this to:

alphac=(exp((Vd/mV-10)/11-(Vd/mV-6.5)/27)/18.975)*(Vd<=50*mV)/ms+2*exp((6.5-Vd/mV)/27)*(Vd>50*mV)/ms : Hz

After making this change it acts exactly how it should.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.