library("bioPN")

onlyR <- F

if (!onlyR) {
  # NOTE: change extension .so to .dll if in Windows
  dyn.load("protocols_propensities.so")
}

source("constants.R")
source("places_transitions.R")
source("model.R")
source("fast_slow.R")
source("initcond.R")
#cat("\nInitial conditions are rounded to zero decimals to start\n")
#cat("as some quantities have to be integers for the stochastic simulation.\n\n")
print(M <- round(M,0))

ect <- 1e-9
delta <- 10*60

runs <- 10

runs.to.plot <- min(30,runs)

if (onlyR) {
  Gy <- 2
  timep <- function() {
    StartTime = round(StartTime,0)
    EndTime = 82*3600
    if (StartTime == 0) {
      y[Gamma] = 0
      y[TNF] = 0
      EndTime = 30*3600
    } else if (StartTime == 30*3600) {
      y[TNF] = 10
      EndTime = 33*3600
    } else if (StartTime == 33*3600) {
      y[Gamma] = Gy
      EndTime = 34*3600
    } else if (StartTime == 34*3600) {
      y[Gamma] = 0
      y[TNF] = 0
    }
    list(EndTime,y)
  }
} else {
  timep <- getNativeSymbolInfo("TNF_before_IR", PACKAGE="protocols_propensities")$address
}

model <- list(pre=pre, post=post, h=h, slow=slow, M=M,
              place=place, transition=transition)
set.seed(19761111)
Sim <- HaseltineRawlings(model, timep, delta, runs, ect)

nLast <- length(Sim$dt)
#ApopSurv <- t(sapply(Sim$run, function(run) {as.matrix(run$M[nLast,c("Aktp","P53pn")])}))
ApopSurv <- t(sapply(Sim$run, function(run) {as.data.frame(run$M)[nLast,c(Aktp,P53pn)]}))

colnames(ApopSurv) <- c("Aktp", "P53pn")
pdf(file="second_example.pdf", width=8, height=7)
par(mar=c(2, 2, 2, 1) + 0.1)
plot(ApopSurv,type="p", pch=".",
     main=paste("Fraction of Apoptotic/Surviving cells (", runs," runs)", sep=""),
     ylim=c(0,8e5),xlim=c(0,5e4))
legend("topleft", "P53pn", bty="n")
legend("bottomright", "Aktp  ", bty="n")
b <- 8e5/5e4
abline(a=0,b=b)
text(1.9e4,4.5e5, paste("Apoptotic:", sum(ApopSurv[,"P53pn"] >= as.numeric(ApopSurv[,"Aktp"]) * b)))
text(3.6e4,4.5e5, paste("Surviving:", sum(ApopSurv[,"P53pn"] < as.numeric(ApopSurv[,"Aktp"]) * b)))

plot(Sim$dt/3600,Sim$run[[1]]$M[[P53pn]],type="l",col="green",ylim=c(0,6e5),
     xlab="",ylab="", main=paste("P53pn, MDM2pn and DSB (first", runs.to.plot, "runs)"),
     panel.first=grid(col="black"))
lines(Sim$dt/3600,Sim$run[[1]]$M[[MDM2pn]],type="l",col="red")
lines(Sim$dt/3600,(Sim$run[[1]]$M[[DSB]])*1000,type="l",col="blue")    
if (runs.to.plot > 1) {
  for (n in 2:runs.to.plot) {
    lines(Sim$dt/3600,Sim$run[[n]]$M[[P53pn]],type="l",col="green")
    lines(Sim$dt/3600,Sim$run[[n]]$M[[MDM2pn]],type="l",col="red")
    lines(Sim$dt/3600,Sim$run[[n]]$M[[DSB]]*1000,type="l",col="blue")
  }
}
old.par <- par(font=2)
legend("topleft",c("P53pn", "MDM2pn", "DSB*1000"),
       col=c("green", "red", "blue"), lty=1)
par(old.par)

plot(Sim$dt/3600,Sim$run[[1]]$M[[AF]],type="l",col="orange",ylim=c(0,9e5),
     xlab="",ylab="", main=paste("AF (first", runs.to.plot, "runs)"),
     panel.first=grid(col="black"))
abline(h=6.5e5, lty=2, col="red")
if (runs.to.plot > 1) {
  for (n in 2:runs.to.plot) {
    lines(Sim$dt/3600,Sim$run[[n]]$M[[AF]],type="l",col="orange")
  }
}

which.places <- c(
                  TNF, Gamma,
                  G_PTEN,
                  G_MDM2,
                  G_P53,
                  G_A20,
                  G_IkBa
                  )
for (plc in which.places) {
  plot.type = "s"
  
  plot(Sim$dt/3600,Sim$run[[1]]$M[[plc]],type=plot.type,main=paste(place[plc], "(first run)"),
       xlab="",ylab="",ylim=c(0, max(as.numeric(lapply(Sim$run, function(run) {max(run$M[[plc]])})))),
       panel.first=grid(col="black"))
}
which.places <- c(
                  PTENt, PTEN,
                  PIP, PIPp,
                  Akt, Aktp,
                  MDM2t, MDM2, MDM2p, MDM2pn,
                  P53t, P53n, P53pn,
                  IKKKn, IKKKa,
                  IKKn, IKKa, IKKi, IKKii,
                  A20t, A20,
                  IkBat, IkBa, IkBap, IkBan, 
                  NFkB, NFkBn,
                  NFkB_IkBa, NFkB_IkBap, NFkB_IkBan,
                  Rec,
                  DSB
                  )
for (plc in which.places) {
  plot.type = "l"
  
  plot(Sim$dt/3600,Sim$run[[1]]$M[[plc]],type=plot.type,main=paste(place[plc],
                                                         "(first", runs.to.plot, "runs)"),
       xlab="",ylab="",ylim=c(0, max(as.numeric (lapply(Sim$run, function(run) {max(run$M[[plc]])})))),
       panel.first=grid(col="black"))
  
  if (runs.to.plot > 1) {
    for (n in 2:runs.to.plot) {
      lines(Sim$dt/3600,Sim$run[[n]]$M[[plc]],type=plot.type)
    }
  }

}    
dev.off()