#############################################################################
# Deterministic model with ecology based on Boots, M. and Y. Haraguchi. 1999. 
# The evolution of costly resistance in host-parasite systems. Am. Nat. 153: 359-370
# Coded by D. George and C. Webb May 28, 2008, updated May 13, 2009
#############################################################################
rm(list=ls())
require(odesolve)

### Parameter values #########################################
pars = c(a= 2, b= 0.012, q= 0.003, g= 1, sigma= 0.01)
tInt = seq(0, 10000, by=1)         # time steps
Initial = c(S= 6, I= 1, B= 0.4)

LinearCostEcol <- function(t, y, parms){
    S <- y[1]
    I <- y[2]
    B <- y[3]
    #Equations for system of ODEs
    Sdot <- (parms[1]*B+parms[2])*S - parms[3]*(S+I)*S - B*S*I
    Idot <- B*S*I - parms[4]*I
    Bdot <- parms[5] * (parms[1]-I)
    OEDs <- c(Sdot, Idot, Bdot)
    list(OEDs)
}

output2 = lsoda(Initial, times= tInt, func= LinearCostEcol, parms= pars)
head(output2)

### Plotting ##########################
par(mfrow=c(3,1))
plot(tInt, output2[,"S"], "l", xlab="", main="Susceptible Population Size", ylab="Susceptible", lwd=1.5, col="dark blue")
plot(tInt, output2[,"I"], "l", xlab="", main="Infected Population Size", ylab="Infected", lwd=1.5, col="dark red")
plot(tInt, output2[,"B"], "l", xlab="Time", main="Evolution of Susceptibility", ylab="Mean Susceptibility (B)", lwd=1.5, col="dark green")