#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

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
m=0.3;      #mutation probability
N=50000;    #Number of events that will occur

#Initial Values  **************************************************
X.initial.size=20;          #Initial susceptible population size
Y.initial.size=3;           #Initial infected population size
W.initial.size=0;           #Initial wait time
S <- rep(0, N)
W <- S
X <- S
Y <- S
X[1] <- X.initial.size
Y[1] <- Y.initial.size
W[1] <- S[1]
loopTime=0

B <- rep(0.3, X.initial.size)  #Initial susceptibility values  

trait=array(0, 1000);  
trait[1:X.initial.size]=B;
traitCumulative=matrix(0,1000,1000)

for( i in 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  ) break
    
    #determine time to next event from exponential distribution (Sojurn time)
    #S[i]=rexp(1,rate=1/sumrates)
    S[i]=rexp(1,rate=sumrates)
    W[i]=W[i-1]+S[i]   #update wait time

    #Decide what event occurred and update populations
    u<-runif(1,min=0,max=1)

    c=length(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=runif(1,min=0,max=1)   #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 while (s1 > cumtotal/Xborn)

       #mutation process
       s2=runif(1,min=0,max=1)
       if(s2 < m){
           B[c+1]=B[selection]+rnorm(1, mean=0, sd=0.01)  #mutation
	     } else{ B[c+1]=B[selection] }    #no mutation
       
       #end susceptible is born
       
       } else{ if(u<(Xborn+Xdie)/sumrates){  #susceptible dies
          X[i]=X[i-1]-1
          Y[i]=Y[i-1]
      
          #decide who dies
          s1=runif(1,min=0,max=1)   #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 while (s1> cumtotal/Xdie) 
            B[selection]=-200
        
    #end susceptible dies
        
       } else{ if(u<(Xborn+Xdie+Yborn)/sumrates){    #infected born
          X[i]=X[i-1]
          Y[i]=Y[i-1]+1
    
    #end infected born
    
       } else{ if(u>(Xborn+Xdie+Yborn)/sumrates){    #infected dies
          X[i]=X[i-1]
          Y[i]=Y[i-1]-1
       } 
  
    #end infected dies

    } #end third else if  "infected dies"
    } #end second else if "infected born"
    } #end first else if "susceptible dies"
    
            
   B=B[B>-1] #shortens B array to the same size as the population (i.e., remove dead individuals)
   #trait[i,]=rep(0,1000)
   #trait[i,1:X[i]]=B; #record trait values
   
   if(length(B)==0) break         
   traitCumulative[i,]=rep(0,1000)
   traitCumulative[i,1:X[i]]=B        
   loopTime=loopTime+1
            
   if(X[i]<=0) break; #end
   if(Y[i]<0)  break; #end                                      
} # end for loop                   

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

par(mfrow=c(2,1))
#figure(1)
plot(c(W,W),c(X,Y),xlab="Time",ylab="Population Size",type="n")
lines(W,Y,type="l",col="dark green",lwd=2)      #infecteds
lines(W,X,type="l",col="blue",lwd=2)            #susceptibles
#legend('Susceptibles','Infecteds')
legend("right", c("Susceptibles","Infecteds"),lwd=c(2,2),col=c("blue","dark green"))

cut=length(W)
numPlots=seq(from=1,to=cut,by=10)
traitCumulative=traitCumulative[1:cut,]
plot(traitCumulative[,cut],W,xlab="Average susceptibility, B",ylab="Time",col="black",cex=1.5,pch=21,bg="grey",xlim=c(0,1))
for(i in 1:length(numPlots)){
  points(traitCumulative[,numPlots[i]],W,col="black",cex=1.5,pch=21,bg="grey")
}  
par(mfrow=c(1,1))