# model of Bahrad Sokhansanj NAR 30, 1817-25 (2002) # see also Free Radical Biology and Medicine 37,422-427 (2004) # Evren Gurkan's matlab code converted into R by Tom Radivoyevitch on 11/4/2005 library(odesolve) e1=48.8; e2=119; e3=19.8; e4=19.8; e5=19.8; e6=39.7; # kcat Km # 1/seconds nM k1=0.00139; K1=1863; # ogg1 excision # values from the 2002 paper # Without Ogg1-Ape Coop. k1=7.55e-3; K1=8.9;# With Ogg1-Ape Coop k2=0.00089; K2=13.7; # ogg1 AP lyase k3=3.08; K3=32.5; # APE endo k4=0.05; K4=130; # APE 3' diesterase k5=0.45; K5=300; # pol beta gap filling k6=0.0022 ; K6=67; # pol delta gap filling k7=0.075*3;K7=500; # With PolBeta-Ape Coop. #k7=0.075;K7=500;# Without PolBeta-Ape Coop. k8=0.39; K8=39; # Fen1 flap endo k9=0.04; K9=100; # lig1 ligation # assume this for a background rate per day Vbkgrnd=20000/(24*3600) # rate created per second from rate per day fbahrad <- function(t, y, p) { v1=y[1]*k1*e1/(y[1]+K1); v2=y[2]*k2; v3=y[2]*k3*e2/(y[2]+K3); v4=y[4]*k4*e2/(y[4]+K4); v5=y[3]*k5*e3/(y[3]+K5); v6=y[6]*k5*e3/(y[6]+K5); v7=y[5]*k7*e3/(y[5]+K7); v8=y[6]*k6*e4/(y[6]+K6); v9=y[3]*k6*e4/(y[6]+K6); v10=y[9]*k8*e5/(y[9]+K8); v11=y[7]*k9*e6/(y[7]+K9); v12=y[8]*k9*e6/(y[8]+K9); v13=y[10]*k9*e6/(y[10]+K9); dy=rep(0,11); dy[1]=-v1+Vbkgrnd; dy[2]=v1-v2-v3; dy[3]=v3-v5-v9; dy[4]=v2-v4; dy[5]=v5-v7; dy[6]=v4-v6-v8; dy[7]=v7-v11; dy[8]=v6-v12; dy[9]=v8+v9-v10; dy[10]=v10-v13; dy[11]=v11+v12+v13; list(dy,c(Va=v12,Vb=v13,Vc=v11))} out1=lsoda(y=rep(0,11),times=seq(-2000,0,1),fbahrad, parms=c(trt=0), rtol=1e-4, atol= rep(1e-4,11)) ny0=out1[nrow(out1),2:12] ny0[1]=200 # assume we add 200 lesions by an acute dose of ionizing radiation out2=lsoda(y=ny0,times=seq(0,3000,1),fbahrad, parms=c(trt=0), rtol=1e-4, atol= rep(1e-4,11)) outs=data.frame(rbind(out1,out2)) outs #outs=data.frame(out1) attach(outs) par(mfrow=c(3,3)) plot(time,Va,type="l",xlab="Time (sec)") plot(time,Vb,type="l",xlab="Time (sec)") plot(time,Vc,type="l",xlab="Time (sec)") # note that long patch handles most of the lesions plot(time,X1,type="l",xlab="Time (sec)") plot(time,X2,type="l",xlab="Time (sec)") plot(time,X3,type="l",xlab="Time (sec)") plot(time,X5,type="l",xlab="Time (sec)") plot(time,X7,type="l",xlab="Time (sec)") plot(time,X11,type="l",xlab="Time (sec)") par(mfrow=c(1,1)) detach(outs)