#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
rm(list=ls())
require(odesolve)

#Models ******************************************************************
linearcostecol=function(t,y,parms){
    #Assign values of y to state variables.
    #This order must be maintained in the inital conditions.
    X1 = y[1]
    X2 = y[2]
    X3 = y[3]

    #Equations for system of ODEs
    x1dot=(a*X3+b)*X1-q*(X1+X2)*X1-X3*X1*X2
    x2dot=X3*X1*X2-g*X2
    x3dot=sigma*(a-X2);

    #The function must return a list containing the derivatives.
    xdot = c(x1dot,x2dot,x3dot)
    return(list(xdot))
}

#parameter values ****************************************************
a= 2
b= 0.012
q= 0.003
g= 1
sigma = 0.01
parms = numeric(5)

#initial conditions ********************************************
x10 = 20
x20 = 1
x30 = 0.4
times=seq(0,2000,by=1) # time steps
x0=c(x10,x20,x30)

output1 = lsoda(x0,times,linearcostecol,parms)

#plotting *******************************************************************
tm=output1[,1]           #time no evolution
popn1=output1[,2]        #population size susceptibles
popn2=output1[,3]		     #population size infected
trait=output1[,4]		     #mean susceptibility

par(mfrow=c(3,1))
plot(tm,popn1,xlab="Time",ylab="Susceptible Population Size",lwd=2,col="red",type="l")
title(main="Susceptible Population Size")
plot(tm,popn2,xlab="Time",ylab="Susceptible Population Size",lwd=2,col="red",type="l")
title(main="Susceptible Population Size")
plot(tm,trait,xlab="Time",ylab="Mean Susceptibility",lwd=2,
  col="red",type="l")
title(main="Evolution of Susceptibility")