# Title:  	CONTINUAL CHANGE IN MATE PREFERENCES
# Source: 	Nature [0028-0836] IWASA yr:1995 vol:377 iss:6548 pg:420 -422 
#---------------------------------------------------------------------------------

# An example of adaptive dynamics with some frequency dependent selection

# The fitness gradients are not given in the paper and have to be computed.

# They are partial derivatives of W_f and W_m towards each of the two individual trait values.
# In the derivative the individual trait value p or t has then to be replaced by the population average trait value.

# ln(W_f) = -b * p^2

# Does not depend on average trait value in the population

# partial derivative for p

# d/dp ln(W_f) = 
# -b * 2  * p

# ln(W_m) = a*p_av*(t - t_av) - c*t^4


# d/dt ln(W_m) =
# a*p_av  - 4*c*t^3

# Depends on average trait value for p (p_av)

# --------------------------------------------------------------------------------

# parameters

a<-0.4 
c<-0.05
b<-.001

# G matrix


G <-diag(0.5,2)

# the covariance is claimed to stabilize
# this result is derived in a different paper

G[1,2]<-G[2,1]<-a*(G[1,1]*G[2,2])/2

G

# mutational bias

u<-c(0,0)




# recurrence of the average trait values #
##########################################

update<-function(z,G,u){
grad1<-
a*z[2]  - 4*c*z[1]^3;

grad2<-
-b * 2  * z[2]

return(z+0.5*G%*%c(grad1,grad2)+u)
}


# simulation #
##############

par(mfrow=c(1,2))

Gen <- 10000
trajectory<-array(0,dim=c(Gen+1,2))

trajectory[1,]<-c(0.5,.05)
for(i in 1:Gen)trajectory[i+1,1:2]<-update(trajectory[i,],G,u)
plot(trajectory,xlab="male trait", ylab="preference")


# include mutation bias
u<-c(-.001,0)

trajectory2<-array(0,dim=c(Gen+1,2))
trajectory2[1,]<-c(0.5,.05)
for(i in 1:Gen)trajectory2[i+1,1:2]<-update(trajectory2[i,],G,u)
plot(trajectory2,xlab="male trait", ylab="preference")

tail(trajectory2)