# As in Fig_6c
# We set the kDsB=0 and kauto=0 at t=480 min (or 8 h)
global 1 t-480.001 {kDSBB=0;kauto=0}

#1.p53
dp53m/dt=sp53+e5*(delay(P53p,tau1)+delay(P53pp,tau1))^n/((delay(P53p,tau1)+delay(P53pp,tau1))^n+Kp53^n)-deltap53*p53m
#2.mdm2
dmdm2m/dt=smdm2-deltamdm2*mdm2m+e1*(delay(P53p,tau1)+delay(P53pp,tau1))^n/((delay(P53p,tau1)+delay(P53pp,tau1))^n+Km^n)
#3.P53
dP53/dt=rp53*delay(p53m,tau2)-uP53*P53-katm1*ATMp*P53+kwip1*P53p*Wip1-kf1*Mdm2*P53+kb1*C1-kf3*Mdmx*P53+kb3*C3-k2*P53*C2-k1*C1
#4.P53p
dP53p/dt=katm1*ATMp*P53-kwip1*P53p*Wip1-kphos*P53p+kdephos*P53pp-uP53p*P53p
#5.Mdm2
dMdm2/dt=rmdm2*delay(mdm2m,tau2)-kd2*Mdm2-katm2*Mdm2*ATMp+kwip2*Mdm2p*Wip1-kf1*Mdm2*P53+kb1*C1-kf2*Mdm2*Mdmx+kb2*C2
#6.C1
dC1/dt=kf1*P53*Mdm2-kb1*C1
#7.wip1
dwip1m/dt=swip1-deltawip1*wip1m+e2*(delay(P53p,tau1)+delay(P53pp,tau1))^n/((delay(P53p,tau1)+delay(P53pp,tau1))^n+Kw^n)
#8.Wip1
dWip1/dt=rwip1*delay(wip1m,tau2)-uWip1*Wip1
#9.ATM* +kauto*ATMp*(1-ATMp)
dATMp/dt=kDSBB*DSB/(DSB+KDSB)*(1-ATMp)-kwip4*Wip1*ATMp-kbasal*ATMp+kauto*(1-ATMp)*ATMp
#10.mdmx
dmdmxm/dt=smdmx-deltamdmx*mdmxm+e3*(delay(P53p,tau1)+delay(P53pp,tau1))^n/(Kx^n+(delay(P53p,tau1)+delay(P53pp,tau1))^n)
#11.Mdmx
dMdmx/dt=rmdmx*delay(mdmxm,tau2)-uMdmx*Mdmx-katm3*Mdmx*ATMp+kwip3*Mdmxp*Wip1-kf2*Mdm2*Mdmx+kb2*C2-kf3*P53*Mdmx+kb3*C3
#12.Mdmxp
dMdmxp/dt=katm3*Mdmx*ATMp-kwip3*Mdmxp*Wip1-k3*Mdmxp*Mdm2p-uMdmxp*Mdmxp
#13.Mdm2p
dMdm2p/dt=katm2*Mdm2*ATMp-kwip2*Mdm2p*Wip1-k4*Mdm2p*Mdm2-uMdm2p*Mdm2p
#14. C2
dC2/dt=kf2*Mdm2*Mdmx-kb2*C2
#15. C3 
dC3/dt=kf3*P53*Mdmx-kb3*C3
#16. P53pp
dP53pp/dt=kphos*P53p-kdephos*P53pp-uP53pp*P53pp
#17. p21
dp21m/dt=e4*delay(P53a,tau1)^n/(delay(P53a,tau1)^n+Kp21^n)-deltap21*p21m
#18. P21
dP21/dt=rp21*delay(p21m,tau2)-uP21*P21-kd21*C2*P21

#19. dP53a
dP53a/dt=Vmax1*P53/(P53+Km1+Km1*Mdm2/Ki1)+Vmax2*P53/(P53+Km2+Km2*Mdmx/Ki2)-uP53a*P53a

# kd2
kd2=uMdm2+DSB/(J+DSB)*kd22
par kd22=0.01, J=0.2

# parameters value
par sp53=5e-4,smdm2=0.002,swip1=0.002,smdmx=0.001
par deltap53=0.03, deltamdm2=0.05, deltawip1=0.05, deltamdmx=0.03, deltap21=0.04
par rp53=0.01, rmdm2=0.04, rwip1=0.02, rmdmx=0.01, rp21=0.02
par uP53=0.03,uMdm2=0.033,uWip1=0.035,uMdmx=0.03,uP21=0.03,uMdm2p=0.1,uMdmxp=0.2,uP53p=0.01, uP53pp=0.004,uP53a=0.001
par e1=0.02, e2=0.014, e3=0.005, e4=0.015,e5=0.02
par k1=0.2, k2=0.01, k3=1.5, k4=0.1, katm1=0.8, katm2=0.02, katm3=0.02, kwip1=1.3, kwip2=0.5, kwip3=0.2, kwip4=1.5
par kDSBB=0.0005, kbasal=0.02, Km=0.16, Kw=0.2, KDSB=200,kauto=0.07
par DSB=3
par kf1=552, kf2=600, kf3=552 , kb1=123.6 , kb2=18 , kb3=123.6
par Kx=1.5
par kphos=0.3, kdephos=0.05
par Kp53=0.017, Kp21=0.017
par n=4, tau1=30, tau2=10, kd21=0.5
par Vmax1=0.0001, Km1=0.000025, Ki1=0.00001, Vmax2=0.001,Km2=0.5, Ki2=0.0001

# initial conditions
init p53m=0.05, mdm2m=0.05, P53=0.0258, P53p=0, Mdm2=0.15, C1=0.06, wip1m=0.04, Wip1=0, ATMp=0, mdmxm=0.01, Mdmx=0.08, Mdmxp=0.01, Mdm2p=0.0178, C2=0.05, C3=0.05, P53a=0.01,p21m=0.05,P21=0.01

# total protein
aux P53t=P53+P53p+P53pp+P53a
aux Mdm2t=Mdm2+Mdm2p
aux Mdmxt=Mdmx+Mdmxp
#@ NPLOT=4,XP=t, YP=P53t, YP2=Mdm2t, YP3=Mdmxt, YP4=p21m, TOTAL=1500, METH=stiff, XLO=0, XHI=1500, YLO=0, YHI=0.4, delay=3060,maxstor=500000,bound=10000
@ XP=t, YP=P53t, TOTAL=1500, METH=stiff, XLO=0, XHI=1500, YLO=0, YHI=0.4, delay=3060,maxstor=500000,bound=10000
done