High frequency stimulation of the Subthalamic Nucleus (Rubin and Terman 2004)

 Download zip file 
Help downloading and running models
Accession:116867
" ... Using a computational model, this paper considers the hypothesis that DBS works by replacing pathologically rhythmic basal ganglia output with tonic, high frequency firing. In our simulations of parkinsonian conditions, rhythmic inhibition from GPi to the thalamus compromises the ability of thalamocortical relay (TC) cells to respond to depolarizing inputs, such as sensorimotor signals. High frequency stimulation of STN regularizes GPi firing, and this restores TC responsiveness, despite the increased frequency and amplitude of GPi inhibition to thalamus that result. We provide a mathematical phase plane analysis of the mechanisms that determine TC relay capabilities in normal, parkinsonian, and DBS states in a reduced model. This analysis highlights the differences in deinactivation of the low-threshold calcium T -current that we observe in TC cells in these different conditions. ..."
Reference:
1 . Rubin JE, Terman D (2004) High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. J Comput Neurosci 16:211-35 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Thalamus; Basal ganglia; Subthalamic Nucleus;
Cell Type(s): Thalamus geniculate nucleus/lateral principal GLU cell; Subthalamus nucleus projection neuron; Globus pallidus neuron;
Channel(s): I Na,t; I T low threshold; I K; I Calcium;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: XPP;
Model Concept(s): Activity Patterns; Parkinson's; Deep brain stimulation;
Implementer(s):
Search NeuronDB for information about:  Thalamus geniculate nucleus/lateral principal GLU cell; I Na,t; I T low threshold; I K; I Calcium;
# rubin_terman_pd.ode :  For Rubin & Terman, JCNS, 2004
# use for normal case (with rubin_terman_pd.ode.setnorm)
# use for parkinsonian case without DBS (with rubin_terman_pd.ode.setpark)
# use rubin_terman_dbs.ode,.ode.set for parkinsonian case with DBS

# stn parameters
p vl=-60,vna=55.,vk=-80.,thetam=30,sm=15
p gl=2.25,gna=37.5,gk=45,tn=1.,th=0.05
p gahp=9.,gca=.5,vca=140.,k1=15.,eps=5e-05
p kca=22.5,thetas=39.,ss=8.,xp=1.,i=0.
p thetah=-39,sh=3.1,thetan=-32.,sn=-8.
p taun0=1,taun1=100.,thn=80.,sigman=26.
p tauh0=1,tauh1=500,thh=57.,sigmah=3.,phi=.75
p thetat=-63.,kt=-7.8,gt=.5,phir=.5
p thetar=-67,kr=2.,taur0=7.1,taur1=17.5,thr=-68,sigmar=2.2
p alpha=5,beta=1.,ab=-30.,gsyn=0.9,vsyn=-100
p rth=.25,rsig=-.07,rho1=.5,a1=.9,hets=0,inp=25
p alphai=1,betai=.05

# gp parameters
p gnag=120.,gkg=30.,gahpg=30.,gtg=.5,gcag=.1,glg=.1
p vnag=55.,vkg=-80.,vcag=120.,vlg=-55.
p thetasg=-57.,ksg=2.,thetas1g=-35.,ks1g=2.
p thetarg=-70.,krg=-2.,taurg=30.
p thetamg=-37.,sigmamg=10.
p thetang=-50.,sigmang=14.
p taun0g=.05,taun1g=.27,thng=-40,sng=-12
p thetahg=-58,sigmahg=-12
p tauh0g=.05,tauh1g=.27,thhg=-40,shg=-12
p k1g=30.,kcag=15.,epsg=0.0001
p phig=1.,deltang=.1,deltahg=.05
p iapp=2,gsyngg=.3,vsyngg=-80.
p gsyng=0.3,vsyng=0,alphag=2,betag=.04,abg=-20
p iappi=5,igl=1.,thresg=0.0,hetg=0.,betagi=.08,kcagi=15
p gsyngi=.3,gsynggi=1.5,alphaggi=1,betaggi=.1

# TC parameters
p vsyntc=-85,gsyntc=.08,asg=200,bsg=.4
p itc=6,shi=-80,shi2=-90,dur=5,dur2=10,period=25
p gnabar=3,gkbar=5,glbar=.05,ena=50,ek=-90,eleak=-70
p gtbar=5,qht=2.5,tadj=1,apr=4,apt=1

# stn functions
sinf(v)=1./(1.+exp(-(v+thetas)/ss))
minf(v)=1./(1.+exp(-(v+thetam)/sm))
hinf(v)=1./(1.+exp((v-thetah)/sh))
ninf(v)=1./(1.+exp((v-thetan)/sn))
taun(v)=taun0+taun1/(1+exp((v+thn)/sigman))
tauh(v)=tauh0+tauh1/(1+exp((v+thh)/sigmah))
rinf(v)=1/(1+exp((v-thetar)/kr))
taur(v)=taur0+taur1/(1+exp((v+thr)/sigmar))
tinf(v)=1/(1+exp((v-thetat)/kt))
rnew(r)=1/(1+exp((r-rth)/rsig))-1/(1+exp(-rth/rsig))

# gp functions
sinfg(vg)=1/(1+exp(-(vg-thetasg)/ksg))
sinf1g(vg)=1/(1+exp(-(vg-thetas1g)/ks1g))
rinfg(vg)=1/(1+exp(-(vg-thetarg)/krg))
minfg(vg)=1./(1.+exp(-(vg-thetamg)/sigmamg))
ninfg(vg)=1./(1.+exp(-(vg-thetang)/sigmang))
taung(vg)=taun0g+taun1g/(1+exp(-(vg-thng)/sng))
hinfg(vg)=1./(1.+exp(-(vg-thetahg)/sigmahg))
tauhg(vg)=tauh0g+tauh1g/(1+exp(-(vg-thhg)/shg))
hv(v)=1./(1.+exp(-v/.001))

# TC functions
inaf(v,m,h)=gnabar*m^3*h*(v-ena)
ikf(v,n)=gkbar*n^4*(v-ek)
ilf(v)=glbar*(v-eleak)
itf(v,mt,ht)=gtbar*mt^2*ht*v
minftc(v)  =  1/(1+exp(-(v+37)/7))
mtinf(v)  = 1/(1+exp(-(v+60)/6.2))
ah(v) =  0.128*exp(-(46+v)/18)
bh(v) =  apr/(1+exp(-(23+v)/5))
tauhtc(v)  =  1/(ah(v)+bh(v))
hinftc(v)  =  1/(1+exp((v+41)/4))
htinf(v)  = 1/(1+exp((v+84)/4))
tauht(v)=(28+apt*exp((v+25)/(-10.5)))

# function for cortical inputs, periodic case:
ftc(t)=itc*hv(sin(6.2831853*(t+shi)/period))*(1-hv(sin(6.2831853*(t+shi+dur)/period)))

# function used to set up windows for calculation of error index:
ftc2(t)=hv(sin(6.2831853*(t+shi2)/period))*(1-hv(sin(6.2831853*(t+shi2+dur2)/period)))

# stn currents
il(v)=gl*(v-vl)
ina(v,h)=gna*(minf(v))^3*h*(v-vna)
ik(v,n)=gk*n^4*(v-vk)
iahp(v,ca)=gahp*(v-vk)*ca/(ca+k1)
ica(v)=gca*((sinf(v))^2)*(v-vca)
it(v,r)=gt*(tinf(v)**3)*(rnew(r)^2)*(v-vca)

isyn1=gsyn*(sg5+sg2)*(v1-vsyn)
isyn2=gsyn*(sg6+sg1)*(v2-vsyn)
isyn3=gsyn*(sg8+sg4)*(v3-vsyn)
isyn4=gsyn*(sg7+sg3)*(v4-vsyn)
isyn5=gsyn*(sg2+sg6)*(v5-vsyn)
isyn6=gsyn*(sg1+sg5)*(v6-vsyn)
isyn7=gsyn*(sg3+sg8)*(v7-vsyn)
isyn8=gsyn*(sg4+sg7)*(v8-vsyn)
isyn9=gsyn*(sg13+sg10)*(v9-vsyn)
isyn10=gsyn*(sg14+sg9)*(v10-vsyn)
isyn11=gsyn*(sg16+sg12)*(v11-vsyn)
isyn12=gsyn*(sg15+sg11)*(v12-vsyn)
isyn13=gsyn*(sg10+sg14)*(v13-vsyn)
isyn14=gsyn*(sg9+sg13)*(v14-vsyn)
isyn15=gsyn*(sg11+sg16)*(v15-vsyn)
isyn16=gsyn*(sg12+sg15)*(v16-vsyn)


# gp currents
itg(vg,rg)=gtg*(sinfg(vg)^3)*rg*(vg-vcag)
inag(vg,hg)=gnag*(minfg(vg)^3)*hg*(vg-vnag)
ikg(vg,ng)=gkg*(ng^4)*(vg-vkg)
iahpg(vg,cag)=gahpg*(vg-vkg)*cag/(cag+k1g)
icag(vg)=gcag*((sinf1g(vg))^2)*(vg-vcag)
ilg(vg)=glg*(vg-vlg)

isyng1=gsyng*(s8+s4+.2*s16)*(vg1-vsyng)
isyng2=gsyng*(s7+s3+.2*s15)*(vg2-vsyng)
isyng3=gsyng*(s1+s5+.2*s9)*(vg3-vsyng)
isyng4=gsyng*(s2+s6+.2*s10)*(vg4-vsyng)
isyng5=gsyng*(s4+s8+.2*s12)*(vg5-vsyng)
isyng6=gsyng*(s3+s7+.2*s11)*(vg6-vsyng)
isyng7=gsyng*(s5+s2+.2*s13)*(vg7-vsyng)
isyng8=gsyng*(s6+s1+.2*s14)*(vg8-vsyng)

isyng9=gsyng*(s16+s12+.2*s8)*(vg9-vsyng)
isyng10=gsyng*(s15+s11+.2*s7)*(vg10-vsyng)
isyng11=gsyng*(s9+s13+.2*s1)*(vg11-vsyng)
isyng12=gsyng*(s10+s14+.2*s2)*(vg12-vsyng)
isyng13=gsyng*(s12+s16+.2*s4)*(vg13-vsyng)
isyng14=gsyng*(s11+s15+.2*s3)*(vg14-vsyng)
isyng15=gsyng*(s13+s10+.2*s5)*(vg15-vsyng)
isyng16=gsyng*(s14+s9+.2*s6)*(vg16-vsyng)


isyngg1=gsyngg*(sg2+sg3)*(vg1-vsyngg)
isyngg2=gsyngg*(sg1+sg5)*(vg2-vsyngg)
isyngg3=gsyngg*(sg4+sg8)*(vg3-vsyngg)
isyngg4=gsyngg*(sg3+sg1)*(vg4-vsyngg)
isyngg5=gsyngg*(sg6+sg7)*(vg5-vsyngg)
isyngg6=gsyngg*(sg5+sg2)*(vg6-vsyngg)
isyngg7=gsyngg*(sg8+sg3)*(vg7-vsyngg)
isyngg8=gsyngg*(sg7+sg4)*(vg8-vsyngg)
isyngg9=gsyngg*(sg10+sg11)*(vg9-vsyngg)
isyngg10=gsyngg*(sg9+sg13)*(vg10-vsyngg)
isyngg11=gsyngg*(sg12+sg16)*(vg11-vsyngg)
isyngg12=gsyngg*(sg11+sg9)*(vg12-vsyngg)
isyngg13=gsyngg*(sg14+sg15)*(vg13-vsyngg)
isyngg14=gsyngg*(sg13+sg10)*(vg14-vsyngg)
isyngg15=gsyngg*(sg16+sg11)*(vg15-vsyngg)
isyngg16=gsyngg*(sg15+sg12)*(vg16-vsyngg)


isyngi[1..16]=gsyngi*si[j]*(vgi[j]-vsyng)

isynggi[1..16]=gsynggi*sggi[j]*(vgi[j]-vsyn)

i V1=-78.,V2=-55.,V3=-78.,V4=-52.
i V5=-78.,V6=-78.,H1=0.2,H2=0.08
i H3=0.2,H4=0.08,H5=0.2,H6=0.2
i N1=0.08,N2=0.4,N3=0.08,N4=0.43
i N5=0.08,N6=0.08,R1=0.67,R2=0.38
i R3=0.67,R4=0.38,R5=0.67,R6=0.67
i CA1=0.44,CA2=0.47,CA3=0.4,CA4=0.4
i CA5=0.44,CA6=0.44,S1=0.001,S2=0.3
i S3=0.001,S4=0.35,S5=0.001,S6=0.31
i VG1=-16.,VG2=-69.,VG3=-6.8,VG4=-69.
i VG5=-8.,VG6=-69.,NG1=0.84,NG2=0.2
i NG3=0.8,NG4=0.2,NG5=0.83,NG6=0.2
i HG1=0.1,HG2=0.7,HG3=0.14,HG4=0.7
i HG5=0.15,HG6=0.7,RG1=0.6,RG2=0.47
i RG3=0.66,RG4=0.47,RG5=0.66,RG6=0.47
i CAG1=0.068,CAG2=0.06,CAG3=0.067,CAG4=0.05
i CAG5=0.067,CAG6=0.05,SG1=0.025,SG2=0.025
i SG3=0.06,SG4=0.025,SG5=0.025,SG6=0.025

i V10=-78.,V9=-512.,V11=-78.,V13=-129.
i V12=-78.,V15=-78.,H10=0.9,H9=0.08
i H11=0.9,H13=0.08,H12=0.9,H15=0.9
i N10=0.08,N9=0.13,N11=0.08,N13=0.1311
i N12=0.08,N15=0.08,R10=0.157,R9=0.118
i R11=0.157,R13=0.118,R12=0.157,R15=0.157
i CA10=0.1313,CA9=0.137,CA11=0.13,CA13=0.13
i CA12=0.1313,CA15=0.1313,S10=0.0010,S9=0.11
i S11=0.0010,S13=0.1112,S12=0.0010,S15=0.1110
i VG10=-1015.,VG9=-159.,VG11=-15.8,VG13=-159.
i VG12=-8.,VG15=-159.,NG10=0.813,NG9=0.9
i NG11=0.8,NG13=0.9,NG12=0.811,NG15=0.9
i HG10=0.10,HG9=0.7,HG11=0.1013,HG13=0.7
i HG12=0.1012,HG15=0.7,RG10=0.15,RG9=0.137
i RG11=0.1515,RG13=0.137,RG12=0.1515,RG15=0.137
i CAG10=0.0158,CAG9=0.015,CAG11=0.0157,CAG13=0.012
i CAG12=0.0157,CAG15=0.012,SG10=0.0912,SG9=0.0912
i SG11=0.015,SG13=0.0912,SG12=0.0912,SG15=0.0912

i vt1=-65.7,ht1=.008577

# stn equations
v[1..16]'=-(il(v[j])+ina(v[j],h[j])+ik(v[j],n[j])+iahp(v[j],ca[j])+ica(v[j])+it(v[j],r[j]))-isyn[j]+hets*[j]+inp
h[1..16]'=phi*( hinf(v[j])-h[j] )/tauh(v[j])
n[1..16]'=phi*( ninf(v[j])-n[j] )/taun(v[j])
r[1..16]'=phir*(rinf(v[j])-r[j])/taur(v[j])
ca[1..16]'=phi*eps*(-ica(v[j])-it(v[j],r[j]) - kca*ca[j])
s[1..16]'=alpha*(1-s[j])*sinf(v[j]+ab)-beta*s[j]
si[1..16]'=alphai*(1-si[j])*sinf(v[j]+ab)-betai*si[j]

# gp equations
vg[1..16]'= -(itg(vg[j],rg[j])+inag(vg[j],hg[j])+ikg(vg[j],ng[j])+iahpg(vg[j],cag[j])+icag(vg[j])+ilg(vg[j]))+iapp-isyngg[j]-isyng[j]+hetg*[j]
ng[1..16]'= deltang*(ninfg(vg[j])-ng[j])/taung(vg[j])
hg[1..16]'= deltahg*(hinfg(vg[j])-hg[j])/tauhg(vg[j])
rg[1..16]'= phig*(rinfg(vg[j])-rg[j])/taurg
cag[1..16]'=epsg*(-icag(vg[j])-itg(vg[j],rg[j]) - kcag*cag[j])
sg[1..16]'=alphag*(1-sg[j])*sinfg(vg[j]+abg)-betag*sg[j]
sggi[1..16]'=alphaggi*(1-sggi[j])*sinfg(vg[j]+abg)-betaggi*sggi[j]

# gpi equations
vgi[1..16]'= -(itg(vgi[j],rgi[j])+inag(vgi[j],hgi[j])+ikg(vgi[j],ngi[j])+iahpg(vgi[j],cagi[j])+icag(vgi[j])+ilg(vgi[j]))+iappi-isyngi[j]-isynggi[j]
ngi[1..16]'= deltang*(ninfg(vgi[j])-ngi[j])/taung(vgi[j])
hgi[1..16]'= deltahg*(hinfg(vgi[j])-hgi[j])/tauhg(vgi[j])
rgi[1..16]'= phig*(rinfg(vgi[j])-rgi[j])/taurg
cagi[1..16]'=epsg*(-icag(vgi[j])-itg(vgi[j],rgi[j]) - kcagi*cagi[j])
sgi[1..16]'=alphag*(1-sgi[j])*sinfg(vgi[j]+abg)-betagi*sgi[j]

# TC equations
vt1'=-(ilf(vt1)+inaf(vt1,minftc(vt1),ht1)+ikf(vt1,.75*(1-ht1))+itf(vt1,mtinf(vt1),htt1))+ftc(t)-gsyntc*sga*(vt1-vsyntc)
htt1'=qht*(htinf(vt1)-htt1)/tauht(vt1)
ht1'=tadj*(hinftc(vt1)-ht1)/tauhtc(vt1)

vt2'=-(ilf(vt2)+inaf(vt2,minftc(vt2),ht2)+ikf(vt2,.75*(1-ht2))+itf(vt2,mtinf(vt2),htt2))+ftc(t)-gsyntc*sgb*(vt2-vsyntc)
htt2'=qht*(htinf(vt2)-htt2)/tauht(vt2)
ht2'=tadj*(hinftc(vt2)-ht2)/tauhtc(vt2)

# additional variables used for error index calculations
coun1'=0
coun2'=0
z[1..2]'=1-hv(z[j]-50)
good[1..2]'=0

p thr0=20
i z1=50,z2=50
i good1=0,good2=0

global 1 {(vt1+thr0)*hv(t-25)*hv(ftc2(t)-0.5)} {coun1=coun1+1}
global 1 {(vt2+thr0)*hv(t-25)*hv(ftc2(t)-0.5)} {coun2=coun2+1}
global 1 {(vt1+10)*hv(t-25)*hv(ftc2(t)-0.5)*hv(z1-40)} {z1=0;good1=good1+1}
global 1 {(vt2+10)*hv(t-25)*hv(ftc2(t)-0.5)*hv(z2-40)} {z2=0;good2=good2+1}

sga=sgi1+sgi2+sgi5+sgi6+sgi9+sgi10+sgi13+sgi14
sgb=sgi3+sgi4+sgi7+sgi8+sgi11+sgi12+sgi15+sgi16

aux sg=ftc(t)
aux sgt=ftc2(t)

aux sa=sga
aux sb=sgb


@ maxstor=100001,dt=.1,total=1000,meth=qualrk,xp=t,yp=v1,xlo=0,xhi=1000,ylo=-80,yhi=20.,bound=50000

done