#Adaptive dynamics simulation from
#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()) #clears values in 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=1000;   #Number of events that will occur

#Initial Values  **************************************************
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 ){

	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]);

    #fitnesses
    fitnessres=r_res-q*(Xeq+Yeq)-B[i-1]*Yeq;
    fitnessmut=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(fitnessmut>fitnessres){
        B[i]=Bmut		#mutant wins
	} else{
        B[i]=B[i-1]} #resident wins
} # end for loop

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

plot(W,B,xlab="Time",ylab="Susceptibility",type="p")