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()