XPP model

This model was converted from XPP ode format to SBML using sbmlutils-0.1.5a6.

# This ODE file reproduces figure 2 from Rempe, Best, Terman J. Math. Biol. 2010
# Notice that as it is written now the simulation contains noise.  If you want 
# to reproduce Figure 2 exactly, turn off noise by making gnoise=0 and gnoiserem=0

# To reproduce figure 3, you will need noise and you'll need to modify the code below 
# to start at different circadian phases.  You'll also probably need MATLAB or some 
# other computing language to process the output of the XPP runs.  

# The default mode is to plot the wake-active, sleep-active, NREM-active, and 
# REM-active populations as in figure 2 panel b.  You can also uncomment the 
# second-to-last line to see the 2-process behavior.  



# Sleep/wake parameters
par eps=3,gamma=5.7,ky=0.01,thetay=0
par epsv=3,gammav=3.77,kv=0.01,thetav=0
par xcirc=-2.1,ghom=5.5
par gsynv=2
par ic0=0,ic1=0
par hmax=1,alphah=18.2,betah=4.2
par gscn=1
par gsyn=5
par const=100


#initialize inp and inpv
par inpbase=3.3,inpvbase=0.45


# ---- REM-NREM stuff ----
par epsr=0.3,gammar=6.2,kr=0.01
par epsn=0.3,gamman=6,kn=0.01
par taur1=1,taur2=2
par taun1=0.5,taun2=1.7
par grem=5,gnrem=0.4
par gevlpo=6.2
par REMscale=11
par inpn=1.9,inpr=1.1
par gA2R=2.5

taur=sinf(xr,kr)*taur2+(1-sinf(xr,kr))*taur1
taun=sinf(xn,kn)*taun2+(1-sinf(xn,kn))*taun1
sinf(v,k)=1/(1+exp(-(v)/k))


#inhibitory currents
inhbr=grem*sn
inhbn=gnrem*sr
ievlpo=gevlpo*xe
iA2R=gA2R*s1

#AMIN-to-REM dynamics
s1'=as1*(1-s1)*sinf(x,0.1)-bs1*s1*(1-sinf(x,0.1))
par as1=2,bs1=0.55

isyn=gsyn*sv
isynv=gsynv*s

iscn=gscn*(xc-xcirc)
par gorx=1
iorx=gorx*iscn*(1-heav(xv))
ihom=ghom*h


tauv=sinf(xv,kv)*tauv2+(1-sinf(xv,kv))*tauv1
tau=sinf(x,ky)*taua2+(1-sinf(x,ky))*taua1

par tauv1=1,tauv2=2
par taua1=1,taua2=2

global 0 t {ic0=gammav-4-inpv}
global 0 t {ic1=gsynv-inpv}


# MonoAminergic group
x'=(3*x-x^3+2-y-isyn+iorx-ihom+inp)*const
y'=eps*(gamma*sinf(x,ky)-y)/tau
s=sinf(x,0.1)


# VLPO group
xv'=(3*xv-xv^3+2-yv-isynv+inpv-iscn+ihom)*const
yv'=epsv*(gammav*sinf(xv,kv)-yv)/tauv
sv=sinf(xv,0.1)


#eVLPO
par ke=0.25,a=2,b=1,c=1.82,kxv=1,kx=0.2,cv=-0.3
xe'=-xe+(c-a*sinf(xv+cv,kxv)-b*sinf(x,kx))
global 0 t {xe=(c-a*sinf(xv,kxv)-b*sinf(x,kx))}


#REM-active population
par rconst=2
xr'=rconst*REMscale*(3*xr-xr^3+3-yr-inhbr+inpr-iA2R+inoiserem)
yr'=REMscale*(epsr*(gammar*sinf(xr,kr)-yr)/taur)
sr=sinf(xr,krn)
par krn=1



#NREM-active population
xn'=rconst*REMscale*(3*xn-xn^3+3-yn-inhbn+inpn-ievlpo+inoisenrem)
yn'=REMscale*(epsn*(gamman*sinf(xn,kn)-yn)/taun)
sn=heav(xn)   
#sn=sinf(xn,0.1)  #either one works


#circadian curve
tc=t-2.5
xc=(0.97*sin((pi/12)*tc)+0.22*sin(2*(pi/12)*tc)+0.07*sin(3*(pi/12)*tc)+0.03*sin(4*(pi/12)*tc)+0.001*sin(5*(pi/12)*tc))


#homeostat
h'=(hmax-h)*(sinf(x,kv))/alphah-h*(1-sinf(x,kv))/betah
#initialize homeostat to what it would normally start at after 
#a normal night of sleep in normal phase
global 0 t {h=0.0978}


#noise stuff
par a1=20,w1=0.01,gnoise=5,gnoiserem=5
ina'=(-ina+a1*normal(0,1))*w1
inv'=(-inv+a1*normal(0,1))*w1
inn'=(-inn+a1*normal(0,1))*w1
inr'=(-inr+a1*normal(0,1))*w1
inoisea=gnoise*ina
inoisev=gnoise*inv
inoiserem=gnoiserem*inr
inoisenrem=gnoiserem*inn



# add some noise to the otherwise constant inputs
inp=inpbase+inoisea
inpv=inpvbase+inoisev



#initial conditions
init x=2,y=6
init xv=-2,yv=0
init xr=-1,yr=-2
init xn=-3,yn=0
init ina=0
init inv=0
init inn=0
init inr=0


#AUXILIARY VARIABLES
# comment out for faster speed
aux ic=iscn+ic1
aux ica=iscn+ic0
aux ih=ihom
#aux noise=inoisea

# first one plots 2-process stuff, next one plots rem, nrem, aminergic
#@ dt=.001,total=48,meth=euler,xp=t,yp=ih,xlo=0,xhi=48,ylo=-3,yhi=3,nplot=3,yp2=ic,yp3=ica,bound=1000000,maxstor=2000001
@ dt=.001,total=48,meth=euler,xp=t,yp=xn,xlo=0,xhi=48,ylo=-3,yhi=3,nplot=3,yp2=xr,yp3=x,bound=1000000,maxstor=2000001,output=output1.dat

done

Terms of use