model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    N <- D + S
    Fr <-  1/(1+N/hR)
    FpS <- 1/(1+N/hS)
    FpD <- 1/(1+N/hD)
    dtS <- (2*Fr-1)*FpS*pS*S - dS*S
    dtD <- 2*(1-Fr)*FpS*pS*S + FpD*pD*D - dD*D
    return(list(c(dtS, dtD)))  
  }) 
}  

# Model by Lander et al.

p <- c(hR=1,hS=1,hD=100,pS=1,pD=0,dS=0.1,dD=0.1)
s <- c(S=1,D=1)
plane(vector=T)
tweak <- "nsol$N=nsol$S+nsol$D"
plane()
plane(xmax=100,ymax=100,vector=TRUE)
run(1000,tweak=tweak)
p["pD"] <- 0
plane(portrait=TRUE)
run()