model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    fI <- pmax(0,1-(I/h)^n)
    dI <- beta*E*I*fI - delta*I
    dE <- b*E*(1 - (E+I)/K) - d*E - beta*E*I*fI
    return(list(c(dI, dE)))  
  }) 
} 
s <- c(I=0.01,E=1)
p <- c(b=1,K=1,d=0.01,beta=2,h=0.5,delta=0.1,n=1)
kappa <- with(as.list(p),K*(1-d/b))
plane()
lines(c(0,kappa),c(kappa,0),lty="dotted")
lines(c(p["h"],p["h"]),c(0,2*kappa),lty="dashed")
#plane(portrait=TRUE,tstep=0.1,grid=5)
run(state=c(I=0.03,E=kappa),tstep=0.1,traject=TRUE)

# Make Figure
# size<-5#inch
# 
# pdf("epithelialA.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")
# plane()
# lines(c(0,kappa),c(kappa,0),lty="dotted")
# lines(c(p["h"],p["h"]),c(0,2*kappa),lty="dashed")
# run(state=c(I=0.03,E=kappa),tstep=0.1,traject=TRUE)
# newton(state=c(I=0,E=kappa),plot=TRUE)
# newton(state=c(I=0.4,E=0.4),plot=TRUE)
# newton(state=c(I=0,E=0),plot=TRUE)
# dev.off()
# 
# model2 <- function(t, state, parms) {
#   with(as.list(c(state,parms)), {
#     fI <- 1/(1+(I/h)^n)
#     dI <- beta*E*I*fI - delta*I
#     dE <- b*E*(1 - (E+I)/K) - d*E - beta*E*I*fI
#     return(list(c(dI, dE)))  
#   }) 
# }  
# 
# s <- c(I=0.01,E=1)
# p <- c(b=1,K=1,d=0.01,beta=2,h=0.2,delta=0.1,n=5);plane(odes=model2)
# lines(c(0,kappa),c(kappa,0),lty="dotted")
# run(odes=model2,traject=TRUE)
# 
# pdf("epithelialB.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")
# plane(odes=model2)
# lines(c(0,kappa),c(kappa,0),lty="dotted")
# run(odes=model2,state=c(I=0.03,E=kappa),traject=TRUE)
# newton(odes=model2,state=c(I=0.5,E=0.5),plot=TRUE)
# newton(odes=model2,state=c(I=0,E=0),plot=TRUE)
# dev.off()
# 
# pdf("epithelialC.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")
# p <- c(b=1,K=1,d=0.01,beta=4,h=0.25,delta=0.3,n=10)
# plane(odes=model2,portrait=TRUE,tstep=0.1,grid=9,legend=FALSE)
# lines(c(0,kappa),c(kappa,0),lty="dotted")
# newton(odes=model2,state=c(I=0.3,E=0.5),plot=TRUE)
# newton(odes=model2,state=c(I=0.2,E=0.1),plot=TRUE)
# newton(odes=model2,state=c(I=0.22,E=0.1),plot=TRUE)
# newton(odes=model2,state=c(I=0,E=0),plot=TRUE)
# dev.off()