model <- function(t, state, parms){
  with(as.list(c(state,parms)),{
    C0 <- state[1]
    Ci <- state[2:(n-1)]
    Cj <- state[1:(n-2)]
    Cn <- state[n]
    R <- RT - sum(state)
    dC0 <- k1*R*L - (k11+k2)*C0
    dCi <-  k2*Cj - (k11+k2)*Ci
    dCn <- k2*state[n-1] - k11*Cn
    return(list(c(dC0,dCi,dCn)))    
  }) 
} 

p <- c(L=1,RT=1,k1=0.1,k11=1,k2=1,theta=0)
n <- 5; s <- rep(0,n); names(s) <- paste("C",seq(0,n-1),sep="") # n = 1 + Nphospho
f <- newton(run(timeplot=FALSE))
continue(f,x="k11",xmin=1e-3,xmax=10,y=n,log="x")
p["L"] <- 10
f <- newton(run(timeplot=FALSE))
continue(f,x="k11",xmin=1e-3,xmax=10,y=n,log="x",add=TRUE)
p["L"] <- 100
f <- newton(run(timeplot=FALSE))
continue(f,x="k11",xmin=1e-3,xmax=10,y=n,log="x",add=TRUE)
lines(c(1e-3,10),c(0.5,0.5),lwd=2,col="blue")

n <- 3; s <- rep(0,n); names(s) <- paste("C",seq(0,n-1),sep="")
p <- c(L=1,RT=1,k1=0.1,k11=1,k2=1)
f <- newton(run(timeplot=FALSE))
continue(f,x="L",xmin=1e-1,xmax=1000,y=n,log="x")
p["k11"] <- 0.5
f <- newton(run(timeplot=FALSE))
continue(f,x="L",xmin=1e-1,xmax=1000,y=n,log="x",add=TRUE)
p["k11"] <- 2
f <- newton(run(timeplot=FALSE))
continue(f,x="L",xmin=1e-1,xmax=1000,y=n,log="x",add=TRUE)
lines(c(1e-3,10),c(0.5,0.5),lwd=2,col="blue")

#----------------------------------------------------------

proof <- function(t, state, parms){
  with(as.list(c(state,parms)),{
    dk11 <- k11
    dL <- L
    dCn <- (RT*L/(k11/k1 + L))*(k2/(k11+k2))^n - theta
    return(list(c(dk11,dL,dCn)))
  })
}

size<-5#inch

pdf("proof.pdf",width=size,height=size)
par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")

colors[3] <- "red"
p <- c(RT=1,k1=0.1,k2=1,n=10,theta=0.0025)
s <- c(k11=1,L=1,Cn=0)
p["n"] <- 0
plane(odes=proof,show="Cn",zero=FALSE,xmin=0.001,xmax=10,ymin=1e-5,ymax=10,log="xy",legend=FALSE)
colors[3] <- "blue"
p["n"] <- 10
plane(odes=proof,show="Cn",zero=FALSE,xmin=0.001,xmax=10,ymin=1e-5,ymax=10,log="xy",legend=FALSE,add=TRUE)
p["n"] <- 100
plane(odes=proof,show="Cn",zero=FALSE,xmin=0.001,xmax=10,ymin=1e-5,ymax=10,log="xy",legend=FALSE,add=TRUE)

dev.off()