Failure of Deep Brain Stimulation in a basal ganglia neuronal network model (Dovzhenok et al. 2013)

 Download zip file 
Help downloading and running models
Accession:148637
"… Recently, a lot of interest has been devoted to desynchronizing delayed feedback deep brain stimulation (DBS). ... This study explores the action of delayed feedback stimulation on partially synchronized oscillatory dynamics, similar to what one observes experimentally in parkinsonian patients. …" Implemented by Andrey Dovzhenok, to whom questions should be addressed.
Reference:
1 . Dovzhenok A, Park C, Worth RM, Rubchinsky LL (2013) Failure of delayed feedback deep brain stimulation for intermittent pathological synchronization in Parkinson's disease. PLoS One 8:e58264 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Basal ganglia;
Cell Type(s): Subthalamus nucleus projection neuron; Globus pallidus neuron;
Channel(s): I Na,t; I T low threshold; I K; I_AHP; I Ca,p;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: XPP; MATLAB;
Model Concept(s): Synchronization; Parkinson's; Deep brain stimulation;
Implementer(s): Dovzhenok, Andrey [andrey.dovzhenok at uc.edu];
Search NeuronDB for information about:  I Na,t; I T low threshold; I K; I_AHP; I Ca,p;
#this is the stimulation protocol with one electrode placed near the 5th STN neuron in the array 
#with a local (w1=0.3,w2=0) or non-local (w1=0.3,w2=0.1) spread of stimulation signal.
#Used to generate Figures 6 and 8 in PLos One paper.

#stimulation parameters
par fr=0.01,Kn=2, t0=500,w1=0.3, Cn=0.04,C1=1

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

# gp parameters
par gnag=120,gkg=30,gahpg=30,gtg=0.5,gcag=0.1,glg=0.1
par vnag=55,vkg=-80,vcag=120,vlg=-55
par thetasg=-57,ksg=2,thetas1g=-35,ks1g=2
par thetarg=-70,krg=-2,taurg=30
par thetamg=-37,sigmamg=10
par thetang=-50,sigmang=14
par taun0g=0.05,taun1g=0.27,thng=-40,sng=-12
par thetahg=-58,sigmahg=-12
par tauh0g=0.05,tauh1g=0.27,thhg=-40,shg=-12
par k1g=30,kcag=3,epsg=0.0055
par phig=1,deltang=0.3,deltahg=0.1
par iapp=6, gsyngg=0,vsyngg=-80
p gsyng=0.4,vsyng=0,alphag=2.,betag=.14,abg=-20
p ka=.33

# stn functions
sinf(v)=1./(1.+exp(-(v+thetas)/ss))
sinfca(v)=1./(1.+exp(-(v+thetas)/ssca))
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))

# 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*((sinfca(v))^2)*(v-vca)
it(v,r)=gt*(tinf(v)**3)*(rnew(r)^2)*(v-vca)

###################################################################
circular structure
###################################################################

isyn1=gsyn*(sg10+sg1+sg2)*(v1-vsyn)
isyn[2..9]=gsyn*(sg[j-1]+sg[j]+sg[j+1])*(v[j]-vsyn)
isyn10=gsyn*(sg9+sg10+sg1)*(v10-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)

###########################################
isyng[1..10]=gsyng*s[j]*(vg[j]-vsyng)

###########################################

#stimulation current

K=Cn*Heav(t-t0)
tau=0.5/fr
tau0=0.5/fr

Istim[1..3]=0.
Istim4=w1*Istim5
Istim5=2.*Kn*(1./(1.+exp(-(K*(delay(x1,tau0+tau)-delay(x1,tau0)))/13.))-0.5)
Istim6=w1*Istim5
Istim[7..10]=0.

aux Stim=Istim5

# LFP of the neuronal population,assumes recording with electrode at 5th STN neuron
Lfp=0.1*(-isyn5+C1*Istim5-w1*(isyn6+isyn4))

aux Muf=Lfp

# stn initial conditions
i V1=-71.2,V2=-66.2,V3=-55.8,V4=-64.6,V5=-61.2,V6=-69.0,V7=-71.4,V8=-69.3,V9=-68.6,V10=-72.0
i H1=0.18,H2=0.35,H3=0.26,H4=0.28,H5=0.19,H6=0.31,H7=0.15,H8=0.23,H9=0.22,H10=0.34
i N1=0.33,N2=0.02,N3=0.1,N4=0.05,N5=0.26,N6=0.05,N7=0.44,N8=0.11,N9=0.11,N10=0.03
i R1=0.51,R2=0.57,R3=0.3,R4=0.63,R5=0.48,R6=0.54,R7=0.42,R8=0.62,R9=0.5,R10=0.61
i CA1=0.53,CA2=0.52,CA3=0.51,CA4=0.51,CA5=0.51,CA6=0.51,CA7=0.5,CA8=0.50,CA9=0.45,CA10=0.45
i S1=0.015,S2=0.,S3=0.,S4=0.,S5=0.003,S6=0.,S7=0.22,S8=0.,S9=0.,S10=0.

# gp initial conditions
i VG1=-66.42,VG2=-59.05,VG3=32.92,VG4=-59.91,VG5=-62.2,VG6=-59.7,VG7=-77.47,VG8=-61.83,VG9=-42.9,VG10=-74.73
i NG1=0.24,NG2=0.34,NG3=0.46,NG4=0.33,NG5=0.28,NG6=0.33,NG7=0.62,NG8=0.29,NG9=0.38,NG10=0.34
i HG1=0.55,HG2=0.52,HG3=0.51,HG4=0.55,HG5=0.58,HG6=0.52,HG7=0.29,HG8=0.59,HG9=0.54,HG10=0.42
i RG1=0.22,RG2=0.1,RG3=0.16,RG4=0.12,RG5=0.19,RG6=0.12,RG7=0.15,RG8=0.18,RG9=0.16,RG10=0.13
i CAG1=0.43,CAG2=0.28,CAG3=0.35,CAG4=0.31,CAG5=0.39,CAG6=0.31,CAG7=0.40,CAG8=0.39,CAG9=0.34,CAG10=0.33
i SG1=0.48,SG2=0.02,SG3=0.33,SG4=0.04,SG5=0.37,SG6=0.04,SG7=0.67,SG8=0.13,SG9=0.29,SG10=0.54

#delay initial conditions
i x1=0.

p Istn=10.,Igpe=0.

aux Mf=0.001*delay(x1,tau0)

om=2.*pi*fr

# harmonic oscillators' equations
x1'=y1
y1'=-1.*om*y1-(om^2)*x1+Lfp

# stn equations
v[1..10]'=-(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]+Istn+C1*Istim[j]
h[1..10]'=phi*( hinf(v[j])-h[j] )/tauh(v[j])
n[1..10]'=phi*( ninf(v[j])-n[j] )/taun(v[j])
r[1..10]'=phir*(rinf(v[j])-r[j])/taur(v[j])
ca[1..10]'=phi*eps*(-ica(v[j])-it(v[j],r[j]) - kca*ca[j])
s[1..10]'=alpha*(1-s[j])*sinf(v[j]+ab)-beta*s[j]

# gp equations
vg[1..10]'= -(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-isyng[j]+Igpe
ng[1..10]'= deltang*(ninfg(vg[j])-ng[j])/taung(vg[j])
hg[1..10]'= deltahg*(hinfg(vg[j])-hg[j])/tauhg(vg[j])
rg[1..10]'= phig*(rinfg(vg[j])-rg[j])/taurg
cag[1..10]'=epsg*(-icag(vg[j])-itg(vg[j],rg[j]) - kcag*cag[j])
sg[1..10]'=alphag*(1-sg[j])*sinfg(vg[j]+abg)-betag*sg[j]


@ dt=.05,total=9999,meth=qualrk,xp=t,yp=v1,xlo=0,xhi=10000,ylo=-80,yhi=20.,bound=500001,maxstor=500001,delay=100.

done

Loading data, please wait...