#Ecological simulation based on
#Boots, M. and Y. Haraguchi. 1999. The evolution of costly resistance in
#host-parasite systems. Am. Nat. 153: 359-370

#Coded by C. Webb and D. George May 28, 2008

rm(list=ls())

#Parameters
q=.003;       #intraspecific competition
g=1;          #background plus disease mortality
a=2;          #shape parameter for cost function
b=0.012;      #shape parameter for cost function
c=0 ;         #shape parameter for cost function
B=0.1         #susceptibility
N=20000;      #Number of events that will occur 

#Initial Values 
X=7;          #Initial susceptible population size
Y=2;          #Initial infected population size
W=0;          #Initial wait time

W.initial.size=0;         #Initial wait time
B.initial.size=0.3;	      #Initial trait value for susceptibility
W <- rep(0, N)
B <- W
S <- W
B[1] <- B.initial.size
W[1] <- S[1]

for( i in 2:N ){
    #Birth and Death rates for X and Y
    Xborn=a*sum(B)+b*X[i-1]-q*(X[i-1]+Y[i-1])*X[i-1]
    Xdie=sum(B)*Y[i-1]
    Yborn=sum(B)*Y[i-1]
    Ydie=g*Y[i-1]
    
    sumrates=Xborn+Xdie+Yborn+Ydie;
    
    if (sumrates <= 0  ) break
    
    #determine time to next event from exponential distribution (Sojurn time)
    S[i]=rexp(1,rate=sumrates)
    W[i]=W[i-1]+S[i]   #update wait time
    
    #Decide what event occurred and update populations
    u<-runif(1,min=0,max=1)
    
    if(u < (Xborn)/sumrates){ #susceptible born
          X[i]=X[i-1]+1
          Y[i]=Y[i-1]
       
    #end susceptible born  
       
    } else{ if(u<(Xborn+Xdie)/sumrates){  #susceptible dies
          X[i]=X[i-1]-1
          Y[i]=Y[i-1]       
       
     #end susceptible dies
  
    } else{ if(u<(Xborn+Xdie+Yborn)/sumrates){    #infected born
          X[i]=X[i-1]
          Y[i]=Y[i-1]+1
    
     #end infected born
    
    } else{ if(u>(Xborn+Xdie+Yborn)/sumrates){    #infected dies
          X[i]=X[i-1]
          Y[i]=Y[i-1]-1
    } 
  
    #end infected dies

    } #end third else if  "infected dies"
    } #end second else if "infected born"
    } #end first else if "susceptible dies"
   
   if(X[i]<=0) break; #end
   if(Y[i]<0)  break; #end 
} # end for loop

#plotting******************************************************
#clip vectors to appropriate sizes... this is a hack...
tm=length(W[W>0])+1
W=W[1:tm]
X=X[1:tm]
Y=Y[1:tm]

par(mfrow=c(1,1))
plot(c(W,W),c(X,Y),xlab="Time",ylab="Population Size",type="n")
lines(W,Y,type="l",col="dark green",lwd=2)      #infecteds
lines(W,X,type="l",col="blue",lwd=2)            #susceptibles
#legend('Susceptibles','Infecteds')
legend("right", c("Susceptibles","Infecteds"),lwd=c(2,2),col=c("blue","dark green"))