model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    B <- S + R + I
    dS <- bB*((1-mu)*S+f*mu*R)*(1 - B/K) - dB*S - beta*S*P
    dR <- bB*(f*(1-mu)*R+mu*S)*(1 - B/K) - dB*R
    dI <- beta*S*P - dI*I
    dP <- b*dI*I - dP*P
    return(list(c(dS, dR, dI, dP)))  
  }) 
}  

p <- c(bB=1,dB=0.01,K=1,beta=1,dI=1,dP=0.1,b=10,f=0.9,mu=0.01)
#with(as.list(c(s,p)),h*(k*N/r-1))
s <- c(S=1,R=0,I=0,P=0.1)
run(100,tweak="nsol$B=nsol$S+nsol$R+nsol$I")