%Adaptive dynamics with ecology 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

close all
clear

%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
m=0.3; %mutation probability

%Initial Values
X=20; %Initial susceptible population size
Y=3;    %Initial infected population size
W=0;   %Initial wait time
B(1:X)=0.3;  %Initial susceptibility values  
trait=0*ones(1,1000);  %value in parentheses is approximately K, carrying capacity
trait(1:X)=B;
N=50000; %Number of events that will occur


for i=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 || sumrates == 0)
        break
    end
    
    %determine time to next event from exponential distribution (Sojurn
    %time)
    S(i)=exprnd(1/sumrates); 
    W(i)=W(i-1)+S(i);   %update wait time
    
    %Decide what event occurred and update populations
    
    u=rand();
    
    [r,c1]=size(B);   %determine size of B
    if(u<Xborn/sumrates) %susceptible born
       X(i)=X(i-1)+1;
       Y(i)=Y(i-1);
       
       %decide who gives birth
       s1=rand();   %gives random number to determine event
       selection=0;
       cumtotal=0; 
       while ( s1 > cumtotal/Xborn)
           selection=selection+1;
           cumtotal=cumtotal+(B(selection)^1.007+a)-q*(X(i-1)+Y(i-1));  
       end
       
       s2=rand();
       %mutation process
       if(s2 < m)
           B(c1+1)=B(selection)+normrnd(0,0.01);  %mutation
       else
           B(c1+1)=B(selection); %no mutation
       end
   
   elseif(u<(Xborn+Xdie)/sumrates)  %susceptible dies
       X(i)=X(i-1)-1;
       Y(i)=Y(i-1);
       
       %decide who dies
       s1=rand();   %gives random number to determine event
       selection=0;
       cumtotal=0; 
       while ( s1 > cumtotal/Xdie)
           selection=selection+1;
           cumtotal=cumtotal+(B(selection)*Y(i-1));
       end
       B(selection)=-200;
      
   elseif(u<(Xborn+Xdie+Yborn)/sumrates)    %infected born
       X(i)=X(i-1);
       Y(i)=Y(i-1)+1;
       
   else                 %infected dies
       X(i)=X(i-1);
       Y(i)=Y(i-1)-1;
   end
   
   B=B(B>-1); %shortens B array to the same size as the population (i.e., remove dead individuals)
   trait(i,:)=ones(1,1000)*0;
   trait(i,1:X(i))=B; %record trait values
   
   if(X(i)<0 || X(i)==0)
       break;
   end
   if(Y(i)<0)
       break;
   end
end

figure(1)
plot(W,X,W,Y)
xlabel('Time')
ylabel('Population Size')
legend('Susceptibles','Infecteds')

figure(2)
plot(trait,W,'.')
ylabel('Time')
xlabel('Average susceptibility, B')
axis([0 1 0 max(W)])

