# model of Hahnfeldt et al (1999) Cancer Research, 59, 4770-4775
# coded by Tom Radivoyevitch on 8/7/2005

library(odesolve)
lambda1=.192;lambda2=0;b=5.85;d=0.00873
y0<-c(V=180,K= 625,g=0)

fphil <- function(t, X, p)
{
names(X)<-c("V","K","g")
lx=as.list(X)
attach(lx)
lp=as.list(p)
attach(lp)

gp = -clr*g
Vp = -lambda1*V*log(V/K)
Kp = -lambda2*K+b*V-d*K*V^.667-e*K*g
XP = c(Vp,Kp,gp)
aux=c(e=e,lambda1=lambda1,b=b,d=d)
detach(lx)
detach(lp)
list(XP,aux) 
}

out1=lsoda(y=y0,times=seq(0,20,1),fphil, parms=c(e=0,clr=0), rtol=1e-4, atol= rep(1e-4,3))
outs=as.data.frame(out1)
attach(outs)
par(mfcol=c(2,2))
plot(time,V,type="l",ylab="Tumor Size (uL)",xlab="Time (days)",main="Control",ylim=c(0,12000))
detach(outs)
days=c(0,4,7,10,13,16,19)
size=c(192,728,2319,3590,5852,8591,10438)
points(days,size)


y0<-c(V=180,K= 625,g=20)
outs=NULL
for (i in 0:13) {
out1=lsoda(y=y0,times=seq(i,i+1,.05),fphil, parms=c(e=0.15,clr=0.38), rtol=1e-4, atol= rep(1e-4,3))
y0=out1[nrow(out1),2:4]+c(0,0,20)
outs=rbind(outs,out1)}
outs=as.data.frame(outs)
attach(outs)
plot(time,V,type="l",ylab="Tumor Size (uL)",xlab="Time (days)",main="Angiostatin (20mg/kg/day)",ylim=c(0,500))
#plot(time,g,type="l",ylab="drug (mg/kg)",xlab="Time (days)",main="Angiostatin (20mg/kg/day)")
detach(outs)
days=c(0,4,7,10,13          )
size=c(176,246,265,136,207  )
points(days,size)


y0<-c(V=180,K= 625,g=30)
outs=NULL
for (i in 0:6) {
out1=lsoda(y=y0,times=seq(2*i,2*(i+1),.05),fphil, parms=c(e=1.3,clr=10.1), rtol=1e-4, atol= rep(1e-4,3))
y0=out1[nrow(out1),2:4]+c(0,0,30)
outs=rbind(outs,out1)}
outs=as.data.frame(outs)
attach(outs)
plot(time,V,type="l",ylab="Tumor Size (uL)",xlab="Time (days)",main="TNP470 (30mg/kg/2day)",ylim=c(0,2500))
#plot(time,g,type="l",ylab="drug (mg/kg)",xlab="Time (days)")
detach(outs)
days=c(0,4,7,10,13 )
size=c(170,381,875,1363,1750 )
points(days,size)


y0<-c(V=180,K= 625,g=20)
outs=NULL
for (i in 0:9) {
out1=lsoda(y=y0,times=seq(i,i+1,.05),fphil, parms=c(e=0.66,clr=1.7), rtol=1e-4, atol= rep(1e-4,3))
y0=out1[nrow(out1),2:4]+c(0,0,20)
outs=rbind(outs,out1)}
outs=as.data.frame(outs)
attach(outs)
plot(time,V,type="l",ylab="Tumor Size (uL)",xlab="Time (days)",main="Endostatin (20mg/kg/day)",ylim=c(0,400))
#plot(time,g,type="l",ylab="drug (mg/kg)",xlab="Time (days)",main="Angiostatin (20mg/kg/day)")
detach(outs)
days=c(0,4,7,10 )
size=c(186,157,135,76 )
points(days,size)