#############################################################################
# Adaptive dynamics 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()) 	# clear memory

### Parameters #####################################
q= 0.003		# intraspecific competition coefficient
g= 1			# background plus disease mortality
a= 2			# shape parameter for cost function
b= 0.012		# shape parameter for cost function
N= 10000		# Number of events that will occur

### Initial Values ############################	   
W <- numeric(N)				# create a vector of zeros of length N
B <- W						# create vector to store susceptibility values as it mutates
B[1] <- 0.4					# Initial trait value for susceptibility

for (i in 2:N){

	W[i]= W[i-1] + 1                                # time vector
	Bmut= B[i-1] + rnorm(1, mean= 0, sd= 0.001)     # find mutant value

    # Linear tradeoff
	r.res <- a*B[i-1] + b
	r.mut <- a*Bmut + b

    # Ecological equilibria evaluated at resident strategy
	Xeq <- g/B[i-1]
	Yeq <- (r.res - q*g/B[i-1]) / (q+B[i-1])

    # Determine fitnesses of Resident & Mutant
	fit.res <- r.res - q*(Xeq+Yeq) - B[i-1]*Yeq
	fit.mut <- r.mut - q*(Xeq+Yeq) - Bmut*Yeq

    # Determining who wins (NOTE: We are assuming that if we have an ESS it is 
    # convergence stable, although we are not testing for this condition. 
    # Technically, you do need to test this additional condition,
    # but the condition holds for the examples we are working with and I
    # wanted to keep the code relatively simple.
    
	if (fit.mut > fit.res) B[i]= Bmut	# Mutant wins
	else B[i]= B[i-1]					# Resident wins
} 	# end for-loop

### Plotting ###############################
### clip vectors to appropriate sizes...
clip= length(W[W > 0]) + 1
W= W[1:clip]
B= B[1:clip]

plot(W, B, "l", xlab= "Time", ylab= "Susceptibility, B", main="Adaptive Dynamics", col= "navy", lwd= "1.5")