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







This file has been produced by sbmlutils.

Terms of use

Copyright © 2017 Matthias Koenig

Redistribution and use of any part of this model, with or without modification, are permitted provided that the following conditions are met:

  1. Redistributions of this SBML file must retain the above copyright notice, this list of conditions and the following disclaimer.
  2. Redistributions in a different form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
This model is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.


Model :

id
name
time
substance
extent
volume
area
length
Access SBML model  L3V1

FunctionDefinitions [5] name math sbo cvterm
max minimum x y x x y y
min maximum x y x x y y
heav heavyside x 0 x 0 0.5 x 0 1 x 0 0
mod modulo x y x y x y x 0 y 0 x y x y
sinf v k 1 1 v k

Parameters [102] name constant value unit derived unit sbo cvterm
eps eps = 3 3.0 None
gamma gamma = 5.7 5.7 None
ky ky = 0.01 0.01 None
thetay thetay = 0 0.0 None
epsv epsv = 3 3.0 None
gammav gammav = 3.77 3.77 None
kv kv = 0.01 0.01 None
thetav thetav = 0 0.0 None
xcirc xcirc = -2.1 -2.1 None
ghom ghom = 5.5 5.5 None
gsynv gsynv = 2 2.0 None
ic0 ic0 = 0 0.0 None
ic1 ic1 = 0 0.0 None
hmax hmax = 1 1.0 None
alphah alphah = 18.2 18.2 None
betah betah = 4.2 4.2 None
gscn gscn = 1 1.0 None
gsyn gsyn = 5 5.0 None
const const = 100 100.0 None
inpbase inpbase = 3.3 3.3 None
inpvbase inpvbase = 0.45 0.45 None
epsr epsr = 0.3 0.3 None
gammar gammar = 6.2 6.2 None
kr kr = 0.01 0.01 None
epsn epsn = 0.3 0.3 None
gamman gamman = 6 6.0 None
kn kn = 0.01 0.01 None
taur1 taur1 = 1 1.0 None
taur2 taur2 = 2 2.0 None
taun1 taun1 = 0.5 0.5 None
taun2 taun2 = 1.7 1.7 None
grem grem = 5 5.0 None
gnrem gnrem = 0.4 0.4 None
gevlpo gevlpo = 6.2 6.2 None
remscale remscale = 11 11.0 None
inpn inpn = 1.9 1.9 None
inpr inpr = 1.1 1.1 None
ga2r ga2r = 2.5 2.5 None
as1 as1 = 2 2.0 None
bs1 bs1 = 0.55 0.55 None
gorx gorx = 1 1.0 None
tauv1 tauv1 = 1 1.0 None
tauv2 tauv2 = 2 2.0 None
taua1 taua1 = 1 1.0 None
taua2 taua2 = 2 2.0 None
ke ke = 0.25 0.25 None
a a = 2 2.0 None
b b = 1 1.0 None
c c = 1.82 1.82 None
kxv kxv = 1 1.0 None
kx kx = 0.2 0.2 None
cv cv = -0.3 -0.3 None
rconst rconst = 2 2.0 None
krn krn = 1 1.0 None
a1 a1 = 20 20.0 None
w1 w1 = 0.01 0.01 None
gnoise gnoise = 5 5.0 None
gnoiserem gnoiserem = 5 5.0 None
x x = 2 2.0 None
y y = 6 6.0 None
xv xv = -2 -2.0 None
yv yv = 0 0.0 None
xr xr = -1 -1.0 None
yr yr = -2 -2.0 None
xn xn = -3 -3.0 None
yn yn = 0 0.0 None
ina ina = 0 0.0 None
inv inv = 0 0.0 None
inn inn = 0 0.0 None
inr inr = 0 0.0 None
s1 0.0 dimensionless None
xe 0.0 dimensionless None
h 0.0 dimensionless None
taur 0.0 dimensionless None
taun 0.0 dimensionless None
inhbr 0.0 dimensionless None
inhbn 0.0 dimensionless None
ievlpo 0.0 dimensionless None
ia2r 0.0 dimensionless None
isyn 0.0 dimensionless None
isynv 0.0 dimensionless None
iscn 0.0 dimensionless None
iorx 0.0 dimensionless None
ihom 0.0 dimensionless None
tauv 0.0 dimensionless None
tau 0.0 dimensionless None
s 0.0 dimensionless None
sv 0.0 dimensionless None
sr 0.0 dimensionless None
sn 0.0 dimensionless None
tc 0.0 dimensionless None
xc 0.0 dimensionless None
inoisea 0.0 dimensionless None
inoisev 0.0 dimensionless None
inoiserem 0.0 dimensionless None
inoisenrem 0.0 dimensionless None
inp 0.0 dimensionless None
inpv 0.0 dimensionless None
ic 0.0 dimensionless None
ica 0.0 dimensionless None
ih 0.0 dimensionless None
t model time 0.0 dimensionless None

Rules [44]   assignment name derived units sbo cvterm
d s1/dt = as1 1 s1 sinf x 0.1 bs1 s1 1 sinf x 0.1 None
d x/dt = 3 x x 3 2 y isyn iorx ihom inp const None
d y/dt = eps gamma sinf x ky y tau None
d xv/dt = 3 xv xv 3 2 yv isynv inpv iscn ihom const None
d yv/dt = epsv gammav sinf xv kv yv tauv None
d xe/dt = xe c a sinf xv cv kxv b sinf x kx None
d xr/dt = rconst remscale 3 xr xr 3 3 yr inhbr inpr ia2r inoiserem None
d yr/dt = remscale epsr gammar sinf xr kr yr taur None
d xn/dt = rconst remscale 3 xn xn 3 3 yn inhbn inpn ievlpo inoisenrem None
d yn/dt = remscale epsn gamman sinf xn kn yn taun None
d h/dt = hmax h sinf x kv alphah h 1 sinf x kv betah None
d ina/dt = ina a1 normal 0 1 w1 None
d inv/dt = inv a1 normal 0 1 w1 None
d inn/dt = inn a1 normal 0 1 w1 None
d inr/dt = inr a1 normal 0 1 w1 None
taur = sinf xr kr taur2 1 sinf xr kr taur1 None
taun = sinf xn kn taun2 1 sinf xn kn taun1 None
inhbr = grem sn None
inhbn = gnrem sr None
ievlpo = gevlpo xe None
ia2r = ga2r s1 None
isyn = gsyn sv None
isynv = gsynv s None
iscn = gscn xc xcirc None
iorx = gorx iscn 1 heav xv None
ihom = ghom h None
tauv = sinf xv kv tauv2 1 sinf xv kv tauv1 None
tau = sinf x ky taua2 1 sinf x ky taua1 None
s = sinf x 0.1 None
sv = sinf xv 0.1 None
sr = sinf xr krn None
sn = heav xn None
tc = t 2.5 None
xc = 0.97 12 tc 0.22 2 12 tc 0.07 3 12 tc 0.03 4 12 tc 0.001 5 12 tc None
inoisea = gnoise ina None
inoisev = gnoise inv None
inoiserem = gnoiserem inr None
inoisenrem = gnoiserem inn None
inp = inpbase inoisea None
inpv = inpvbase inoisev None
ic = iscn ic1 None
ica = iscn ic0 None
ih = ihom None
t = time None

Events [4] name trigger priority delay assignments sbo cvterm
e0 t 0
initialValue = False
persistent = True
ic0 = gammav 4 inpv
e1 t 0
initialValue = False
persistent = True
ic1 = gsynv inpv
e2 t 0
initialValue = False
persistent = True
xe = c a sinf xv kxv b sinf x kx
e3 t 0
initialValue = False
persistent = True
h = 0.0978