%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

close all
clear

%Parameters
%q intraspecific competition coefficient
%g background plus disease mortality
%a shape parameter for cost function
%b shape parameter for cost function
%c shape parameter for cost function

%linear tradeoff
q=0.003; g=1; a=2; b=0.012; c=0;

%Initial Values
B=0.3;           %Initial susceptibility value
N=1000;        %Number of timesteps
W=1;            %Initial time
for i=2:N
    W(i)=i; %time vector
    Bmut=B(i-1)+normrnd(0,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;
    else
        B(i)=B(i-1);
    end
 
end

figure(1)
plot(B,W,'.')
ylabel('Time')
xlabel('Susceptibility, B')

