%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 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
B=0.1   %susceptibility

%Initial Values 
X=7; %Initial susceptible population size
Y=2;    %Initial infected population size
W=0;   %Initial wait time

N=20000; %Number of events that will occur


for i=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 || 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();
    
   if(u<Xborn/sumrates) %susceptible born
       X(i)=X(i-1)+1;
       Y(i)=Y(i-1);
       
   elseif(u<(Xborn+Xdie)/sumrates)  %susceptible dies
       X(i)=X(i-1)-1;
       Y(i)=Y(i-1);
       
   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
   
   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')
