#############################################################################
# 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 D. George and C. Webb May 28, 2008, updated May 13, 2009
#############################################################################
rm(list=ls())

#Parameters
q= 0.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 (used in AdaptiveDynamicsEcology)
B= 0.5			# 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 <- numeric(N)				# create a vector of zeros of length N
#B <- W						# used in AdaptiveDynamicsEcology
S <- W
W[1] <- S[1]

for(i in 2:N) {
    # Birth and Death rates for X and Y
    Xborn= (a*B+b)*X[i-1] - q*(X[i-1] + Y[i-1])*X[i-1]
    Xdie= B*X[i-1]*Y[i-1]
    Yborn= B*X[i-1]*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...
clip= length(W[W > 0]) + 1
W= W[1:clip]
X= X[1:clip]
Y= Y[1:clip]

plot(c(W,W), c(X,Y), xlab="Time", ylab="Population Size", type="n") # Setup empty plot
lines(W, X, type="l", col="navy",lwd= 1.5)            # susceptibles
lines(W, Y, type="l", col="dark red",lwd= 1.5)        # infecteds
legend("topleft", c("Susceptibles", "Infecteds"), bg='gray88', lwd=c(2,2), col= c("navy","dark red"))