I'm trying to write a script to demonstrate the effect of selection in a population. The problem that I have is that it is not realistic in the sense that not only the mean would change for directional selection, but the Skewness of the distribution would change as well. How would that be possible to code that and show the movement of the population going toward the maximum fitness peak?

# Define the normal distribution
slect <- function(x,mu,sigma,constant) {
  temp = exp(constant*((x-mu)/sigma)^2)

# Multiply 2 functions 
# Set some parameters 
ylim = c(0,2)
# X limits 
from = -4
to = 7
# Sequence of X 
x = seq(from,to,length.out=100)
# Means 
muinit = 1
# Number of generations 
n.gen = 10
# Variance 
sigma.0 = 1.1
sigma = sigma.0

# Initialize population 
init0 = function(x)slect(x,mu = muinit,sigma = sigma.0,constant = -.5)
# How the population will move because of selection 
popmov = function(x)slect(x,mu = mu1,sigma = sigma,constant = -.5)
# Final fitness function 
finl = function(x)slect(x,mu = mu2,sigma = sigma,constant = -.5)

# Empty data 
data.all = NULL
gen = 1
for(gen in 1:n.gen){
  # system("sleep .2")
  # png(paste0("Desktop/fittest/test",gen,".png"))
  int1 = curve(expr = init0,from = from, to = to, col = "red",
               ylab = "Fitness", xlab = "Phenotype", ylim = ylim, main = paste("Fitness result gen:",gen))

  if (gen ==1) {
    # Result of the first generation 
    h0 = Multiply(init0,finl)
    h1 = curve(expr = h0,from = from, to = to, col = "black", lty = 3, add = TRUE)
    finl1 = curve(expr = finl,from = from, to = to, col = "black", lty = 2,lwd = 2, add = TRUE)
    popmov1 = curve(expr = popmov,from = from, to = to, col = "red",lwd = 2,add = TRUE)
    h = h0
  if (gen>1) {
    mu1 = mean(c(mu1,mu2))
    popmov1 = curve(expr = popmov,from = from, to = to, col = "red",lwd = 2,add = TRUE)
    h = Multiply(h,finl)
    h1 = curve(expr = h,from = from, to = to, col = "black", lty = 3, add = TRUE)
    finl1 = curve(expr = finl,from = from, to = to, col = "black", lty = 2,lwd = 2, add = TRUE)
  popmov1$type = "popmovement"
  finl1$type = "final"
  h1$type = "movement"
  int1$type = "initial"
  popmov1$it = gen
  int1$it = gen
  finl1$it = gen
  h1$it = gen
  df1 = as.data.frame.list(int1)
  df1.1 = as.data.frame.list(popmov1)
  df2 = as.data.frame.list(finl1)
  df3 = as.data.frame.list(h1)
  df4 = rbind(df1,df1.1,df2,df3)
  data.all =c(data.all,list(df4))
  # dev.off()

enter image description here

  • 1
    $\begingroup$ Someone here may be able to answer your question, but if not I would try SE Bioinformatics, which might have been a better bet in the first place. $\endgroup$ – David May 20 at 16:31
  • $\begingroup$ Modeling for a single genotype parameter or a single phenotype parameter is (IMHO) very unrealistic. I know that doesn't address your question, but it may be relevant to your goal. $\endgroup$ – S. McGrew May 20 at 21:11
  • $\begingroup$ @S.McGrew You are right! But my goal here is to have a "simple" case to be visualized (adding traits would start making the thing in 3d and further make it even harder to visualize). But as a first approximation, I wonder how to "show" the effect of a fitness landscape on the population. Therefore, I'm pretty sure that not only the mean is changing, but also the skew of the distribution. This is the core of my question. $\endgroup$ – M. Beausoleil May 20 at 23:41
  • 2
    $\begingroup$ If your selective "force" works equally on all parts of the population, you will keep the shape of the original distribution. In real life, though, the "leading edge" of the population will procreate faster than the "trailing edge". Account for that and I'm pretty sure you will get the skew you're looking for. $\endgroup$ – S. McGrew May 21 at 1:44
  • $\begingroup$ @S.McGrew, This is it, but how to 'account for that' is my question! I don't know how to simulate this part. Now I'm only modifying the mean of the distribution (the ` mu1=mean(c(mu1,mu2))` in the script). But there should be something modifying the shape and skew of the distribution (and Kurtosis if the actual density of the population changes) as well. $\endgroup$ – M. Beausoleil May 21 at 12:24

There are at least several different approaches that can work. Without examining your code, it seems that the population is described as a distribution of phenotypes. You need a way to allow that distribution to change in a fairly complicated way. So if you define a function to describe the change something like this, you should get results closer to what you're expecting:

  1. define distribution as N(p), where (p) is phenotypic distance from the original peak of the distribution, and N is the number of population members at that distance p.
  2. Define D(pf-p) as the distance from the new fitness landscape peak to a point in the distribution that change.
  3. Define a forcing function that depends on D. This function will either decrease or increase the number of population members at each point along the phenotype distribution, depending on D. If you want to get fancy, it can also depend on the number of population members in the previous iteration at that point in the distribution.
  4. If you want the population size to stay constant, you should make sure the forcing function adds and subtracts the same number of members from the population.
  5. I would suggest using a forcing function that gets stronger as the phenotype gets closer to the new fitness landscape peak, and make sure the phenotype distribution is nowhere zero when you begin.

The above function should in effect move population members toward the new fitness peak, though may tend to sharpen the population peak. To fix that you can introduce another function that "blurs" the distribution, like randomly moving members to nearby phenotype values.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.