# model of Ciliberto, A., Novak, B. and Tyson, J. J. Steady states and oscillations in the p53/Mdm2 network. Cell Cycle 4, 488-93 (2005).

# XPP ODE file (obtained from Andrea Ciliberto) converted into R (by hand) by Tom Radivoyevitch on 8/6/2005

library(odesolve)
ks2p=0.0015;ks2=0.006;kd2p=0.01;kd2pp=0.01;
kph=0.05;kdeph=6;ki=14;ko=0.5;
ks53=0.055;kd53=8;kd53p=0.0055;kf=8.8
kr=2.5;kDNA=0.18;kdDNA=0.017

m=3;Js=1.2;Jdam=.2;Vratio=15;
J=0.01

Jdna=1
ampl=1 


fciliberto <- function(t, X, p)
{
names(X)<-c("p53t","p531","p532","m2n","m2c","m2Pc","dam")
lx=as.list(X)
attach(lx)
if (p["trt"]==1) IR=ampl  else IR=0

p53=p53t-(p531+p532)
m2t=m2c+(1/Vratio)*m2n+m2Pc
kd2=kd2p+kd2pp*dam/(dam+Jdam)

p53tp=ks53-kd53p*p53t-kd53*p532
p531p=kf*m2n*p53+kr*p532-p531*(kr+kf*m2n)-kd53p*p531
p532p=kf*m2n*p531-p532*kr-p532*(kd53p+kd53)
m2np=Vratio*(ki*m2Pc-ko*m2n)-kd2*m2n
m2cp=ks2p + ks2*p53t^m/(Js^m+p53t^m) - kd2p*m2c + kdeph*m2Pc - (kph*m2c/(J+p53t))
m2Pcp=kph*m2c/(J+p53t) - kdeph*m2Pc - ki*m2Pc + ko*m2n - kd2p*m2Pc
damp=kDNA*IR - kdDNA*p53t*dam/(dam+Jdna)
XP = c(p53tp,p531p, p532p, m2np, m2cp,m2Pcp,damp );
detach(lx)
list(XP,c(p53=p53,m2t=m2t))}

y0<-c(p53t=0.07,p531=.02,p532=.01,m2n=.33, m2c=.12,m2Pc=.01,dam=0)
out1=lsoda(y=y0,times=seq(-20,0,1),fciliberto, parms=c(trt=0), rtol=1e-4, atol= rep(1e-4,7))
ny0=out1[nrow(out1),2:8]
out2=lsoda(y=ny0,times=seq(0,10,1),fciliberto, parms=c(trt=1), rtol=1e-4, atol= rep(1e-4,7))
ny0=out2[nrow(out2),2:8]
out3=lsoda(y=ny0,times=seq(10,1200,1),fciliberto, parms=c(trt=0), rtol=1e-4, atol= rep(1e-4,7))
outs=data.frame(rbind(out1,out2,out3))
attach(outs)

par(mfcol=c(4,2))
plot(time,p53t,type="l",xlab="Time (min)")
plot(time,p531,type="l",xlab="Time (min)")
plot(time,p532,type="l",xlab="Time (min)")
plot(time,m2n,type="l",xlab="Time (min)")
plot(time,m2c,type="l",xlab="Time (min)")
plot(time,m2Pc,type="l",xlab="Time (min)")
plot(time,dam,type="l",ylab="DNA damage",xlab="Time (min)")
par(mfrow=c(1,1))

detach(outs)