model <- function(t, state, parms) { with(as.list(c(state,parms)), { E <- 1/(1 + (B/h)^n) dB <- s0 + s1*E - d*B return(list(dB)) }) } solveH <- function(H, parms) { with(as.list(parms), { B <- 3e9/d E <- 1/(1 + (B/H)^n) dB <- s0 + s1*E - d*B return(dB) }) } p <- c(s0=1e9,s1=9e9,d=1/120,h=2.5e11,n=5,m=1) s <- c(B=0) run() H <- uniroot.all(solveH,c(1e11,1e12),parms=p) print(c(SolvedH=H)) p["h"] <- H newton(run()) with(as.list(p),curve(1/(1+(B/h)^n),from=0,to=5e11,xname="B",ylim=c(0,1),ylab="EPO",col="red",lwd=2)) with(as.list(p),curve(s0+s1*E^m,from=0,to=1,xname="E",ylim=c(0,1e10),ylab="RBC Production",col="red",lwd=2)) with(as.list(p),curve({E<-1/(1+(B/h)^n);s0+s1*E^m},from=0,to=5e11,xname="B",ylim=c(0,1e10),ylab="Source",col="red",lwd=2)) # Aging: decrease s1 n-fold p young <- newton(run()) p["s0"] <- p["s0"]/2 p["s1"] <- p["s1"]/2 old <- newton(run()) old/young p["s0"] <- p["s0"]*2 p["s1"] <- p["s1"]*2 # make figures size <- 5 #inch pdf("epoA.pdf",width=size,height=size) par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F) with(as.list(p),curve(1/(1+(B/h)^n),from=0,to=5e11,xname="B",ylim=c(0,1),ylab="EPO",col="red",lwd=2)) lines(c(3e9*120,3e9*120),c(0,1)) dev.off() pdf("epoB.pdf",width=size,height=size) par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F) with(as.list(p),curve({E<-1/(1+(B/h)^n);s0+s1*E^m},from=0,to=5e11,xname="B",ylim=c(0,1e10),ylab="Source",col="red",lwd=2)) lines(c(3e9*120,3e9*120),c(0,1e10)) dev.off()