Roff D.A. Modeling Evolution: An Introduction to Numerical Methods
Подождите немного. Документ загружается.
> mean(Output[500:MaxGen,2])
[1] 1.017443
> sd(Output[500:MaxGen,2])
[1] 0.04910655
The trait value evolves toward the value predicted by the optimality approach (1)
and clearly fluctuates about this value (mean 2SD encloses 1). Obviously the
optimality approach is better for efficiently determining the equilibrium value,
but gives no idea of how long this might take or how much the trait value might
fluctuate. The quantitative genetic approach can address these questions and
herein lies its value in the present case.
4.5 Scenario 3: Directional selection using an IVC model
In this scenario we analyze directional selection on a single trait. For the general
quantitative trait the breeder’s equation is adequate if the population size is large.
However, for small populations there will be stochastic variation which will
reduce the rate of change. Formulae to calculate the expected variance about
the response to selection have been worked out for the typical cases (see Roff
[1997], pp. 137–143) but an individual-based quantitative genetic model can be
useful in the case of unusual modes of selection or with traits, such as threshold
traits, that are not normally distributed on the expressed scale. The model
0
0.9 1.0 1.1
Trait value
1.2 1.3 1.4 1.5
500
Generation
1,000 1,500 2,000
Figure 4.3 Time trace for Scenario 2.
248 MODELING EVOLUTION
examined here assumes selection on a threshold trait, but it can be easily modified
to accommodate other types of traits or selection scenarios.
4.5.1 General assumptions
1. The trait under selection has two phenotypic expressions but is determined by
an underlying normally distributed liability.
2. Only one morph is allowed to contribute to the next generation.
4.5.2 Mathematical assumptions
1. Males and females do not differ in the expression of the trait. They will,
nevertheless, be considered separately.
2. Mating between selected individuals is at random.
3. Selected individuals do not differ in their fitness.
4. The two morphs are coded as 0 and 1, with the morph coded as 1 being that
which is selected. This coding means that the morph code also corresponds to
the morph fitness and the mean genetic value of the parents of the next
generation is given by equation (4.18), where fitness, W, is either 0 or 1.
4.5.3 Analysis
In the present coding I have assumed that the trait is wing morph with the short-
winged morph being favored (i.e., short-winged individuals are coded as 1 and
long-winged individuals are coded as 0). The function SELECTION calculates the
new genetic mean liability value using equation (4.18). Although not required, the
phenotypic values of the liability are also passed to the selection function: this is
done here so that the present coding can be more easily modified for selection on a
normally distributed trait that is phenotypically expressed.
At the start of the main program the initial proportion of the two morphs is set
and the required threshold value, Z.LW, calculated assuming, without loss of
generality, a normal distribution with a mean of 0 and a standard deviation of 1
(use the R function qnorm). Male and female short-winged proportions are shown
on the same graph with the values labeled as M or F.
R CODE:
rm(list¼ls()) # Remove all objects from memory
SELECTION <- function(Male, Female) # Function to calculate new
mean value
{
# Mean genetic value of males, mu.M, and females, mu.F
# Col1 ¼ Phenotypic value, col 2 ¼ Genetic value, col 3 ¼ Morph code
mu.M <- sum(Male[,3]*Male[,2])/sum(Male[,3])
mu.F <- sum(Female[,3]*Female[,2])/sum(Female[,3])
mu <- (mu.M þ mu.F)/2 # New population mean
return(mu)
}
GENETIC MODELS 249
######################## Main program ########################
set.seed(100) # Initialize random number generator
N <- 100 # Set population size
Prop.LW <- 0.85 # Set initial proportion LW
Z.LW <- -qnorm(Prop.LW) # Threshold. Values greater are LW
MaxGen <- 10 # Number of generations
Output <- matrix(0,MaxGen,4) # Create file for output
h2 <- 0.5 # Set heritability of liability
Vp <- 1 # Set Phenotypic variance
Va <- Vp*h2 # Calculate Additive genetic variance
Ve <- Vp-Va # Calculate Environmental variance
mu <- 0 # Initial mean genetic liability
SDa <- sqrt(Va) # SD of Va
SDe <- sqrt(Ve) # SD of Ve
for (Igen in 1:MaxGen) # Iterate over generations
{
# Generate Genetic and environmental values using normal distribu-
tion
GM <- rnorm(N, mean¼mu, sd¼SDa) # Genetic values of males
GF <- rnorm(N, mean¼mu, sd¼SDa) # Genetic values of females
EM <- rnorm(N, mean¼0, sd¼SDe) # Environmental values of
males
EF <- rnorm(N, mean¼0, sd¼SDe) # Environmental value of females
PM <-GMþ EM # Phenotypic value of males
PF <-GFþ EF # Phenotypic value of females
# Calculate wing morphs by comparing liability to threshold
Male.Morph <- matrix(1,N) # Set all initially to SW (¼1)
Male.Morph[PM > Z.LW] <- 0 # Set LW to 0
Female.Morph <- matrix(1,N) # Set all initially to SW (¼1)
Female.Morph[PF > Z.LW] <
- 0 # Set LW to 0
#
Combine phenotypic and genetic values
Male <- cbind(PM, GM, Male.Morph )
Female <- cbind(PF, GF, Female.Morph)
# Store data
Output[Igen,1] <- Igen # Generation
Output[Igen,2] <- mean(PMþPF)/2 # Mean liability
Output[Igen,3] <- sum(Male.Morph)/N # Proportion of SW males
Output[Igen,4] <- sum(Female.Morph)/N # Proportion of SW females
# Calculate new mean genetic value by applying fitness criterion
mu <- SELECTION(Male,Female)
} # End of Igen loop
par(mfrow¼c(2,2)) # Divide graphics page into quadrats
# Plot proportion of SW males, and SW females over generation on same
graph
250 MODELING EVOLUTION
plot(Output[,1], Output[,3],pch¼“M”, xlab¼‘Generation’,
ylab¼‘Proportion Short Wings’)
lines(Output[,1], Output[,3]) # Males
points(Output[,1],Output[,4], pch¼“F”) # Females
lines(Output[,1], Output[,4]) # Females
plot(Output[,1], Output[,2], xlab¼ ‘Generation’, ylab¼‘Mean
liability’)
lines(Output[,1], Output[,2])
OUTPUT: (Figure 4.4)
Note that the response to selection is not constant at either the expressed level or
on the underlying liability scale. The reason for this is that the selection intensity
decreases as the proportion of the selected morph increases in the population. For
experimental data showing this process see Roff (1990).
4.6 Scenario 4: Directional selection using an IL model
The previous scenario is here repeated using an individual locus model. The basic
approach is the same as outlined in Section 1.3. This model is a simplified version
of that described in Roff (1994b, 1998a, b).
4.6.1 General assumptions
1. The trait under selection has two phenotypic expressions but is determined by
an underlying normally distributed liability.
2. Only one morph is allowed to contribute to the next generation.
Generation
246810
Generation
Proportion Short Wings
Mean liability
2
0.2
–2.5 –2.0 –1.5 –1.0 –0.5 0.0
0.4 0.6 0.8 1.0
F
46810
M
M
M
M
M
M
M
M
M
M
F
F
F
F
F
F
F
F
F
Figure 4.4 Results for Scenario 3. Response by females labeled “F” and that of males by “M.”
GENETIC MODELS 251
4.6.2 Mathematical assumptions
1. Males and females do not differ in the expression of the trait. They are not
separately distinguished in this model.
2. Mating between selected individuals is at random and can include selfing.
3. Selected individuals do not differ in their fitness.
4. The two morphs are coded as 0 and 1, with the morph coded as 1 being that
which is selected. This coding means that the morph code also corresponds to
the morph fitness and the mean genetic value of the parents of the next
generation is given by equation (4.18), where fitness, W, is either 0 or 1.
5. The phenotypic value is the sum of allelic values and 100 loci, which take
values of 0 or 1, plus an environmental deviation.
4.6.3 Analysis
In Section 1.3 a bivariate model was described: the coding below is a simplified
version of this model. The initial population is started with a random selection of
alleles but the threshold for morph expression is set to 0.85 for long-winged
individuals, as in the previous scenario. Prior to the selection the population is
run for 10 generations to establish an equilibrium distribution of genotypes. The
switch to a selection regime is determined using an “if” statement. A more
efficient method is shown in Scenario 7.
R CODE:
rm(list¼ls()) # Clear memory
SELECTION <- function(Morph, Genotype)
{
Selected <- Genotype[Morph¼¼1,] # 1¼SW, 0¼LW
return(Selected)
} # End function
################## Function Mutation ##################
MUTATION <- function(X, N.mutations, Total.loci, N.ind)
{
# Apply mutation by randomly selecting N.mutations
Row <- ceiling(runif(N.mutations, min¼0, max¼Total.loci))
Temp <- matrix(X)
Temp[Row] <- (abs(Temp[Row]-1))
X <- matrix(Temp,N.ind,N.loci)
return(X)
} # End function
################## Function Gamete ##################
GAMETE <- function(X, G.loci, N.loci) # Pick loci for gamete
{
Y <- sample(x¼X, size¼G.loci, replace¼FALSE) # Random G.loci
from N.loci
return(Y)
}
252 MODELING EVOLUTION
##################### Main Program ######################
set.seed(100) # Initialize random number generator
P <- 0.0001 # Probability of mutation at a locus
N.Pop <- 1000 # Population size at each generation
G.loci <- 100 # Loci per gamete
N.loci <- G.loci*2 # Loci per individual
Total.loci <- N.Pop*N.loci # Total number of loci in population
H2 <- 0.5 # Heritability
Vg <- 0.5*G.loci # Additive genetic variance
Ve <- (1-H2)*Vg/H2 # Environmental variance
SD.E <- sqrt(Ve) # Environmental SD
Prop.LW <- 0.85 # Set initial proportion LW
# Set Threshold value. Values greater than Z.LW are LW
Z.LW <- qnorm(1-Prop.LW, mean¼G.loci, sd¼sqrt(VgþVe))
# Generate matrix of individuals in which
# rows hold individuals while columns hold loci. Allelic values are
1 and 0
# Randomly generate Total.loci number of loci with values of 0 & 1
Dl <- round(runif(Total.loci))
# Genetic composition of individuals
Genotype <- matrix(Dl, N.Pop, N.loci)
Maxgen <- 30 # Number of generations simulation runs
Output <- matrix(0,Maxgen,5) # Allocate space for output
Output[,1] <- seq(from¼1, to¼Maxgen) # First col¼generation
for (Igen in 1:Maxgen) # Iterate over generations
{
Env.X <- rnorm(N.Pop, mean¼0, sd¼SD.E) # Environmental de-
viations
# Phenotypic values of liability
Phenotype <- rowSums(Genotype) þ Env.X
Vg <- var(rowSums(Genotype)) # Calculate the genetic
variance
Vp <- var(Phenotype) # Phenotypic variance
H2 <- Vg/Vp # heritability
Morph <- matrix(1,N.Pop) # Set morphs initially to SW
Morph[Phenotype > Z.LW] <
- 0 # Change relevant individuals
to
LW
Prop.SW <- sum(Morph)/N.Pop # Proportion SW
Output[Igen,2] <- Prop.SW # Store proportion SW
Output[Igen,3] <- H2 # Store heritability
#################### Apply Selection ####################
if (Igen < 10) # No selection until after generation 10
{
Parents <- Genotype # No selection
}else
GENETIC MODELS 253
{
Parents <- SELECTION(Morph, Genotype) # Apply selection
}
N.Parents <- nrow(Parents) # Number of parents
# Form next Generation
# Apply Mutation
# Mean number of mutations in population
lambda <- P*N.Parents*N.loci
# Number of mutations using a Poisson distribution
N.mutations <- rpois(1,lambda)
Output[Igen,4] <- N.mutations # Store number of mutations this
generation
# Apply function MUTATION to generate mutant loci
Genotype <- MUTATION(Genotype, N.mutations, Total.loci, N.Pop)
# Mating
# Produce gametes for female offspring
# Select from each row G.loci at random
# We do not distinguish individual loci
# Note that this creates a matrix of G.loci rows and N.Females
columns
# The matrix is transposed to produce the required matrix
Gametes <- apply(Parents, 1, GAMETE, G.loci, N.loci)
Gametes <- t(Gametes) # Convert to proper matrix
# Produce N.Pop offspring by selecting at random with replacement
Output [Igen,5]<- N.Parents
# Get N.Pop gametes from “females”
n <- seq(1, N.Parents)
G.Index <- sample(x¼n, size¼ N.Pop, replace¼TRUE)
F.Gametes <- Gametes[G.Index,]
# Get N.Pop gametes from “males”
G.Index <- sample(x¼n, size¼ N.Pop, replace¼TRUE)
M.Gametes <- Gametes[G.Index,]
# New Genotypes
Genotype <- cbind(F.Gametes, M.Gametes) # Combine gametes
}
par(mfrow¼c(2,2)) # Divide graph page into four quadrats
plot(Output[,1], Output[,2], xlab¼‘Generation’, ylab¼‘Propor-
tion SW’, type¼‘l’)
plot(Output[,1], Output[,3], xlab¼‘Generation’, ylab¼‘Herit-
ability’,type¼‘l’)
plot(Output[,1], Output[,4], xlab¼‘Generation’, ylab¼‘Nos of
mutations’,type
¼‘l’)
plot(Output[,1],
Output[,5], xlab¼‘Generation’, ylab¼‘Nos of
Parents’,type¼‘l’)
OUTPUT: (Figure 4.5)
254 MODELING EVOLUTION
The general response is the same as in Scenario 3. The primary difference is
that selection leads to a loss in genetic variation, though the effect is slight
and would be very hard to experimentally detect. Also note that after a brief
period following the initiation of selection the population returns very close
to its starting value.
4.7 Scenario 5: A quantitative genetic analysis of the Ricker
model
This scenario investigates the influence of density-dependent selection in the
evolution of a parameter of the Ricker function. It is an individual-based quantita-
tive genetic model of Scenario 3 of Chapter 3.
4.7.1 General assumptions
1. The organism is semelparous.
2. Recruitment is governed by a density-dependence function that allows for
cyclical or chaotic population dynamics.
Generation
Nos of mutations
Proportion SW
Nos of Parents
Heritability
0
5
200 400 600 800 1,000
0.42
0.2 0.4 0.6 0.8 1.0
0.46 0.50
10 15 20 25
5 10152025
30
Generation
0 5 10 15 20 25
30
Generation
0 5 10 15 20 25
30
Generation
0 5 10 15 20 25
30
Figure 4.5 Results for Scenario 4.
GENETIC MODELS 255
3. The parameters of the recruitment function are related such that the density-
independent component is negatively related to the density-dependent compo-
nent.
4. These parameters are inherited.
4.7.2 Mathematical assumptions
1. Population at time t þ 1 is a Ricker function of the population at time t:
N
tþ1
¼ N
t
ae
bN
t
ð4:34Þ
2. The parameter a is a measure of density-independent recruitment whereas b is
a measure of the density-dependent effect: increases in a increase recruitment
but increases in b decrease recruitment by increasing the density-dependent
component, e
bN
t
. Thus a positive functional relationship between a and b is
indicative of a trade-off between the two recruitment components. For this
scenario I shall assume the relationship
b ¼ 0:001a ð4:35Þ
3. The parameter a is a quantitative genetic trait. By the previous assumption, the
correlation between a and b is 1. Note that the genetic correlation could still be
one but that the phenotypic value might still be influenced by a separate
environmental value, in which case we would have
b ¼ 0:001a þ e ð4:36Þ
where e is a random normal variable with mean zero. For simplicity, and to
maintain comparison with the clonal model of Chapter 3, I shall set e ¼ 0.
4. Fitness is the number of offspring left by each individual. Thus individual
fitness, W
i
, which also equates to cohort size, N
i
,is
W
i
¼ N
i
¼ a
i
e
b
i
N
Pop
ð4:37Þ
and the new mean trait value and population size are (subscripts for generation
omitted for simplicity)
m
a
¼
X
N
Pop
i¼1
W
i
G
a;i
X
N
Pop
i¼1
W
i
N
Pop
¼
X
N
Pop
i¼1
W
i
ð4:38Þ
where N
Pop
is the total population size and G
a,i
is the genetic value of the ith
individual.
256 MODELING EVOLUTION
4.7.3 Analysis
The selection function SELECTION requires both the X matrix (consisting of the
phenotypic and genetic values) and the total population size, N. Population size
must be an integer and so is rounded to the nearest integer value. The simulation
is run for 10,000 generations and the mean trait value is calculated after the first
2,000 generations, which appears to be the approximate number of generations
required to achieve the equilibrium value, given the starting value, which is quite
far removed from the optimal value.
R CODE:
rm(list¼ls()) # Remove all objects from memory
SELECTION <- function(X,N) # Function to calculate new mean value
{
BETA <- X[,1]*0.001 # Beta
X.Fitness <- X[,1]*exp(-BETA*N) # Fitness
mu <- sum(X.Fitness*X[,2])/sum(X.Fitness) # New mu
N <- round(sum(X.Fitness)) # Popn size
return(c(mu, N)) # Return values
} # End of selection function
##################### Main program #####################
set.seed(100) # Initialize random number generator
N <- 100 # Set population size
MaxGen <- 10000 # Number of generations
Output <- matrix(0,MaxGen,3) # Create file for output
h2 <- 0.5 # Set heritability
Vp <- .1 # Set Phenotypic variance
Va <- Vp*h2 # Calculate Additive genetic variance
Ve <- Vp-Va # Calculate Environmental variance
mu <- 10 # Trait mean genetic value
SDa <- sqrt(Va) # SD of Va
SDe <- sqrt(Ve) # SD of Ve
for (Igen in 1:MaxGen) # Iterate over generations
{
# Generate Genetic and enviro nmenta l values using normal distribution
GX <- rnorm(N, mean¼mu, sd¼SDa) # Genetic values
EX <- rnorm(N, mean¼0, sd¼SDe) # Environmental values
PX <-GXþ EX # Phenotypic values
# Combine phenotypic and genetic values
X <- cbind(PX,GX)
Output[Igen,1] <- Igen # Store generation
Output[Igen,2] <- mean(PX) # Store mean phenotype
Output[Igen,3] <- N # Store popn size
GENETIC MODELS 257