Roff D.A. Modeling Evolution: An Introduction to Numerical Methods
Подождите немного. Документ загружается.
Trait.X <- 2*(Pxy*X.loci þ Pc*C.loci) # Mean value of X
Trait.Y <- 2*(Pxy*Y.loci þ S*Pc*C.loci) # Mean value of Y
# Genetic variance of iX/i and Y
VarGX <- 2*(X.loci*Pxy*(1-Pxy)þ C.loci*Pc*(1-Pc))
VarGY <- 2*(Y.loci*Pxy*(1-Pxy)þ C.loci*Pc*(1-Pc))
CovGXY <- 2*C.loci*Pc*(1-Pc) # Genetic covariance
Rg <- S*CovGXY/sqrt(VarGX*VarGY) # Genetic correlation
print(c(Rg, VarGX, VarGY, Trait.X, Trait.Y)) # Print out values
# Calculate the environmental correlation
Re <- (Rp - Rg*sqrt(h2.X*h2.Y))/sqrt((1-h2.X)*(1-h2.Y))
# Check that this Re is possible
if (abs(Re)>1) stop (c(“Re not possible”))
# Environmental Variances and Standard deviations
Ve.X <- (1-h2.X)*VarGX/h2.X # Environmental variance for X
SDe.X <- sqrt(Ve.X) # Environmental SD for X
Ve.Y < (1-h2.Y)*VarGY/h2.Y # Environmental variance for Y
SDe.Y <- sqrt(Ve.Y) # Environmental SD for Y
CovE <- Re*SDe.X*SDe.Y # Environmental covariance
Ematrix <- matrix(c (V e.X ,Cov E, Co vE,V e. Y) ,2,2) # Covariance
matrix
# Nos of loci in each category
Nx.Alleles <- ASSIGN.LOCI(X.loci,N.Pop, Pxy)# Alleles uniqueto X
Ny.Alleles<- ASSIGN.LOCI(Y.loci, N.Pop, Pxy) # Alleles uniqueto Y
Nc.Alleles <- ASSIGN.LOCI(C.loci, N.Pop, Pc) # Alleles commontoX&Y
# Now make three matrices for loci in individuals
G.Xmatrix <- matrix(Nx.Alleles, N.Pop, 2*X.loci) # X composition
G.Ymatrix <- matrix(Ny.Alleles, N.Pop, 2*Y.loci) # Y composition
G.Cmatrix <- matrix(Nc.Alleles, N.Pop, 2*C.loci) # C composition
The above coding generates the three matrices. A sample output for G.Xmatrix is
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]
1 1 0 1 1 1100 1
[2,] 1 1 1 1 1 0111 1
[3,] 1 1 1 1 1 1111 1
[4,] 1 1 1 1 1 1111 1
[5,] 0 1 0 1 1 1101 1
[6,] 1 1 1 1 1 1111 1
Note that there are 6 rows, corresponding to 6 individuals, and 10 columns,
corresponding
to 5 diploid loci. We can get the actual genetic variances using
the two R functions rowSums, which calculates the sum of the values in each row
(which is the genetic value of the individual), and var, which calculates the
variances among these sums:
238 MODELING EVOLUTION
# Get actual genotypic values
G.X <- rowSums(G.Xmatrix) þ rowSums(G.Cmatrix) # X Genotypic
values
VarGX <- var(G.X) # Vg for X
G.Y <- rowSums(G.Ymatrix) þrowSums(G.Cmatrix) # Y Genotypic values
VarGY <- var(G.Y) # Vg for Y
print(c(VarGX, VarGY))
For a population size of 1,000 the output from the above code is [1] 2.724881
2.594794, as compared to the expected 2.7426 for both traits. To create the
phenotypic values we add an environmental value to each individual. The
environmental values are correlated and hence we use the function mvrnorm to
generate an N.Pop x 2 matrix of environmental values the first column giving the X
valuesandthesecond columngivingthe Yvalues (remember to load the library MASS):
# Create phenotypic values
Env <- mvrnorm(n¼N.Pop, mu¼c(0,0), Sigma¼Ematrix) # Environmen-
tal values
P.X <- G.X þ Env[,1] # Vector of X phenotypes
P.Y <- G.Y þ Env[,2] # Vector of Y phenotypes
VarPX <- var(P.X) # Phenotypic variance of X
VarPY <- var(P.Y) # Phenotypic variance of Y
The actual heritabilities are thus given by
h2.X <- VarGX /VarPX # Heritability of X
h2.Y <- VarGY/VarPY # Heritability of Y
print(c(h2.X, h2.Y)) # Print results
For the same population of 1,000, as used above, the output is [1] 0.4708531
0.2582277, the expected values being 0.5 and 0.25. The correlations are readily
obtained using the R function cor:
Rg <- cor(G.X, G.Y) # Genetic correlation
Rp <- cor(P.X, P.Y) # Phenotypic correlation
print(c(Rg, Rp))
which gives [1] 0.4947939 0.1973116, where the expected values are 0.5 and
0.25, respectively.
The next process to consider is that of selection. Here I shall assume a threshold
regime in which individuals with an X trait value greater than T0 are selected. In the
coding below T0 is set at the initial mean trait value (Trait.X):
SELECTION <- function(Phenotype, Genotype, T0) # Selection
function
{
Selected <-
Genotype[Phenotype>T0,] # Selection
return(Selected) # Return matrix of selected genotypes
} # End function
GENETIC MODELS 239
In the main program the above function is called three times to select from the three
genotypic matrices according to the values of the phenotypes, in this case the X
phenotype:
ParentX <- SELECTION(P.X, G.Xmatrix, Trait.X) # Select from the X matrix
ParentY <- SELECTION(P.X, G.Ymatrix, Trait.X) # Select from the Y matrix
ParentC <- SELECTION(P.X, G.Cmatrix, Trait.X) # Select from the C matrix
ParentX # Print selecte d matrix
ParentY # Print selecte d matrix
ParentC # Print selecte d matrix
In theaboveexampleselection involves only a singletrait: selectionon multiple traits
can also be done by, for example, using sequential culling or an index. Starting with a
population size of 6 (as in the coding above) the following individuals were selected:
> ParentX # Print selected matrix
[,1]
[,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 0 1 1 1100 1
[2,] 1 1 1 1 1 0111 1
[3,] 1 1 1 1 1 1111 1
[4,] 1 1 1 1 1 1111 1
[5,] 0 1 0 1 1 1101 1
[6,] 1 1 1 1 1 1111 1
> ParentY# Print selected matrix
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 1 1 1 1 1111 1
[2,] 0 1 1 1 0 1111 1
[3,] 1 0 0 1 1 1111 1
[4,] 1 1 1 1 1 1011 1
[5,] 1 0 1 1 1 1111 1
[6,] 1 1 1 1 1 1011 1
> ParentC # Print selected matrix
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 1 1 0 0 1
[2,] 0 0 0 0 1 0
[3,] 0 1 0 1 0 0
[4,] 0 1 0 0 0 0
[5,] 0 1 1 0 0 0
[6,] 0 0 1 1 1 0
The third process to consider is that of offspring production. Given that each
locus
is independent and of equal value (0 or 1) I shall simplify the problem of
forming a gamete pool by ignoring the distinction among loci (i.e., that the first
G.loci represent one set and the second half of the row represents the other set)
and drawing from each individual G.loci from each category at random. For
example, suppose the set of loci for an individual with 5 X loci is given by the
vector
240 MODELING EVOLUTION
Column number 1 2 3 45678910
100001100 0
To exactly mimic sexual reproduction we would select at random between
alleles at positions 1 and 6, 2 and 7, 3 and 8, 4 and 9, 5 and 10. A simpler
approach that should be adequate unless loci differ in their contributions is to
select at random, without replacement, 5 alleles from the set and assign them to
the gamete pool:
GAMETE <- function(X, G.loci) # Pick loci for gamete pool
{
Y <- sample(x¼X, size¼G.loci, replace¼FALSE) # Randomly select G.loci
return(Y)
}
This function is called from the main program using the R function apply. Note
that the resulting matrix has to be transposed to conform to the correct matrix
structure:
# Form Gamete pool
GameteX <- apply(ParentX, 1, GAMETE, X.loci)
GameteX <- t(GameteX) # Convert to proper matrix
GameteY <- apply(ParentY, 1, GAMETE, Y.loci)
GameteY <- t(GameteY) # Convert to proper matrix
GameteC <- apply(ParentC, 1, GAMETE, C.loci)
GameteC <- t(GameteC) # Convert to proper matrix
From the gamete pool we now select 2*N.Pop gametes which are combined at
random to form the next generation of individuals. To do this we first create a
sequence from 1 to N.Parents, where N.Parents is the number of selected
individuals. We next sample at random, with replacement, from this sequence
2*N.Pop values, which become the index values defining the gametes to be
combined to give the next generation. This matrix is called S.GameteJ, where
J ¼ X, Y,orC and consists of 2*N.Pop rows and columns equal to the number of
loci in the relevant category. Individuals are formed by combining the first N.Pop
rows with the second set of rows.
N.Parents <- nrow(ParentX) # Number of available parents
n <- seq(1, N.Parents) # sequence 1 to N.Parents
# Get 2*N.Pop random indices with replacement
G.Index <- sample(x¼n, size¼2*N.Pop, replace¼TRUE)
# Get gametes from gamete pool
S.GameteX <- GameteX[G.Index,]
S.GameteY <- GameteY[G.Index,]
S.GameteC <- GameteC[G.Index,]
# Next generation
GENETIC MODELS 241
n1 <- N.Popþ1
n2 <- 2*N.Pop
G.Xmatrix <- cbind(S.GameteX[1:N.Pop,], S.GameteX[n1:n2,])
G.Ymatrix <- cbind(S.GameteY[1:N.Pop,], S.GameteY[n1:n2,])
G.Cmatrix <- cbind(S.GameteC[1:N.Pop,], S.GameteC[n1:n2,])
At this point we introduce mutation into the scenario. One method is to iterate
over all loci and use a random number generator to determine if the allele at the
locus changes. Such an approach is very time consuming. An alternate approach is
to first generate the number of mutations occurring in this particular set of loci.
Given that the probability of a mutation at a locus is very small the total number of
mutations will follow a Poisson distribution:
PðxÞ¼
e
lx
l!
ð4:29Þ
where x is the number of mutations and l is the mean, given as l ¼ P
M
N
Pop
2n
loci
,in
which P
M
is the per locus mutation rate, N
Pop
is the population size, and n
loci
is the
number of loci in the relevant category. The distribution of the number of
mutations is given by the R function rpois and thus for each generation the
number of mutations can be generated by the call
N.mutations <- rpois(1,lambda)
To select and convert alleles we pass the matrix G.Xmatrix, G.Ymatrix,or
G.Cmatrix to the use function MUTATION, which does the following:
1. Convert the matrix of alleles, called X in the function, into a column vector,
Temp, of length 2*N. inds*G. loci¼T. loci, where G. loci is X. loci, Y.
loci, or C. loci and N. inds is the number of individuals (here N. Pop):
Temp <- matrix(X) # Convert X to a vector
2. Generate N.mutations random integers between 1 and T. loci using the
uniform generator R function runif and the R function ceiling:
Row <- ceiling (runif (N.mutations, min¼0, max¼T. loci))
3. This set of integers represents the rows of the vectors that are affected by
mutation. In principle it is possible to pick the same row twice, though this is
highly unlikely given the low probability of mutation (say 10
4
or less). At each
of these rows the value is changed, 0 to 1 or 1 to 0:
Temp[Row] <- (abs(Temp[Row]-1))
4. The vector is now converted back into a matrix and passed back to the main
program:
X <- matrix(Temp, N.inds, 2*G.loci)
242 MODELING EVOLUTION
The complete function (where Pmut is the mutation probability) is
MUTATION <- function(X, Pmut, G.loci, N.inds)
{
T.loci <- N.inds*2*G.loci # Total number of alleles
lambda <- Pmut*T.loci # Mean number of mutations in population
N.mutations <- rpois(1,lambda) # Number of mutations
Row <- ceiling(runif(N.mutations, min¼0, max¼T.loci))
Temp <- matrix(X) # Convert matrix to a vector
Temp[Row] <- (abs(Temp[Row]-1)) # Change relevant row entries
X <- matrix(Temp, N.inds, 2*G.loci) # Convert back to a matrix
return(X)
} # End function
This completes a single generation of selection. Multiple generations are
simulated by the appropriate inclusion of an iteration loop. The complete pro-
gram is discussed in scenarios 4 and 7.
4.2 Summary of scenarios
Scenario 1: Illustration of the use of a PVC model to plot the trajectory of two
traits under stabilizing selection
Scenario 2: Usage of an IVC model to solve the life history model examined in
Scenario 1 of Chapter 2
Scenario 3: Demonstration of the use of an IVC model to analyze the effect of
directional selection on a single trait
Scenario 4: Illustration of the use of an individual locus model to analyze the
effect of directional selection on a single trait
Scenario 5: An IVC model of Scenario 3 of Chapter 3, namely density-dependent
selection in the evolution of a parameter in the Ricker function
Scenario 6: Analysis of directional selection on a trade-off using an IVC model
Scenario 7: Analysis of directional selection on a trade-off using an individual
locus model
4.3 Scenario 1: Stabilizing selection on two traits using a PVC
model
This scenario explores the consequences of stabilizing selection acting on two
traits as expressed by equation (4.15).
GENETIC MODELS 243
4.3.1 General assumptions
1. Both traits are quantitative with non-zero heritabilities.
2. Selection is stabilizing.
3. Generations do not overlap.
4.3.2 Mathematical assumptions
1. The change in trait values is governed by equation (4.15).
2. For the purposes of this example the heritabilities are set at 0.2 and 0.4, the
phenotypic correlation at 0.4, the genetic correlation at 0.15, and the pheno-
typic variances at 1 and 2. The diagonal elements of the W matrix are set at
2 and 3 and the off-diagonal elements at 0, giving no correlational selection.
The optimum trait values are both set at 2 and the initial trait values at 10.
4.3.3 Analysis
R CODE:
rm(list¼ls()) # Clear workspace
h2 <- c(0.2,0.4) # Set heritabilities
Vp <- c(1,2) # Set phenotypic variances
Rp <- 0.4 # Set phenotypic correlation
Ra <- 0.15 # Set genetic correlation
# Check that Re is possible
Re <- (Rp-Ra*sqrt(h2[1]*h2[2]))/sqrt((1-h2[1])*(1-h2[2]))
if(abs(Re) >1)
{
print (c(“problem with Re”, Re))
stop
}
Va <- h2*Vp # Using h2 ¼Va/Vp
Covp <- Rp*sqrt(Vp[1]*Vp[2]) # Using r ¼ Cov/SD1SD2
Cova <- Ra*sqrt(Va[1]*Va[2]) # Using r ¼ Cov/SD1SD2
Gmatrix <- matrix( c(Va[1],Cova,Cova,Va[2]),2,2) # G matrix
Pmatrix <- matrix(c(Vp[1],Covp,Covp,Vp[2]),2,2) # P matrix
Theta <- c(2,2) # Optimum trait values
Maxgen <- 100 # Number of generations
par(mfrow¼c(2,2)) # Divide graphic page
Wmatrix <- matrix(c(2,0,0,3), 2,2) # Set the W matrix
Trait <- matrix(0,Maxgen,2) # Pre-assign space for trait values
Trait[1,] <- 10 # Initial trait values
for (Igen in 2:Maxgen) # Iterate over generations
{
# Delta z
Delta.Z <- Gmatrix%*%solve(WmatrixþPmatrix)%*%(Trait[Igen-
1,]-Theta)
244 MODELING EVOLUTION
Trait[Igen,] <- Trait[Igen-1,]- Delta.Z # New trait value
} # End of Igen loop
# Set axis values for graphing
min.y <- min.x <- min(Trait); max.y <- max.x <- max(Trait)
# Plot by generation
Generation <- seq(from¼1, to¼Maxgen)
plot(Generation, Trait[,1], ylim <- c(min.y, max.y), xlim¼c(0,
Maxgen), ylab¼‘Trait’, type¼‘l’)
lines(seq(from¼1,to¼Maxgen),Trait[,2],lty¼2)
# Plot Trait 2 on Trait 1
plot(Trait[,1], Trait[,2])
lines(Trait[,1],Trait[,2])
print( c(Rp,Ra,Re))
eigen(Wmatrix)$values
OUTPUT: (Figure 4.2)
> print( c(Rp,Ra,Re))
[1] 0.400000 0.150000 0.516113
> eigen(Wmatrix)$values
[1] 3 2
Both traits eventually reach their respective optima, although it takes quite a few
generations.
4.4 Scenario 2: Stabilizing selection using an IVC model
This scenario is an individual-based quantitative genetic model of Scenario 1 of
Chapter 2.
Generation
0 20406080100
2
24
Trait [ , 2]
Trait [ , 1]
6810
24
Trait
6810
46810
Figure 4.2 Output from Scenario 1 showing the time course of stabilizing selection on two
traits (left) and their joint evolution (right).
GENETIC MODELS 245
4.4.1 General assumptions
1. The organism is semelparous.
2. Fecundity, F, increases with body size, x.
3. Survival, S, decreases with body size, x.
4. Fitness, W, is a function of fecundity and survival.
5. Body size is inherited as a quantitative trait.
4.4.2 Mathematical assumptions
1. Fecundity increases linearly with body size:
F ¼ a
F
þ b
F
x ð4:30Þ
where a
F
and b
F
are constants.
2. Survival decreases linearly with body size:
S ¼ a
S
b
S
x ð4:31Þ
3. Fitness, W, is the expected lifetime reproductive success, R
0
, given as the
product of fecundity and survival:
W ¼ R
0
¼ FS
¼ða
F
þ b
F
xÞða
S
b
S
xÞ
¼ a
F
a
S
b
F
b
S
x
2
þða
S
b
F
a
F
b
S
Þx
ð4:32Þ
The above equation describes a parabola that is concave down. As in the optimality
analysis parameter values are set at a
F
¼ 0, b
F
¼ 4, a
S
¼ 1, and b
S
¼ 0.5. The fitness
equation can now be written as
W ¼ 0 1 4 0:5x
2
þð1 4 0 0:5Þx
¼2x
2
þ 4x
ð4:33Þ
which has a maximum at x ¼ 1.
4.4.3 Analysis
Because it is not possible for individuals to have negative body sizes it may be
necessary to either eliminate these from the population or set them at some mini-
mal value. In the present example this possibility is highly unlikely as the optimal
body size is 10 phenotypic standard deviations from zero. However, to be sure, I set
all negative fecundities, which would correspond to negative body sizes, to zero,
giving such individuals zero fitness. I also insert a check on survival to ensure that it
is always positive. The calculation of the new genetic mean is done within the
function SELECTION. The matrix X passed to SELECTION contains the phenotypic
values in the first column and the genetic values in the second. The long-term mean
phenotype is calculated ignoring the first 500 generations of selection.
R CODE:
rm(list¼ls()) # Remove all objects from memory
246 MODELING EVOLUTION
SELECTION <- function(X) # Function to calculate new mean value
{
As <-1;Bs<- 0.5; Af <-0;Bf<- 4 # Parameter values
Survival <- As-Bs*X[,1] # Survival
Survival[Survival<0] <- 0 # Check on sign
Fecundity <-AfþBf*X[,1] # Fecundity
Fecundity[Fecundity<0] <- 0 # Check on sign
X.Fitness <- Survival*Fecundity # Fitness
mu <- sum(X.Fitness*X[,2])/sum(X.Fitness)
return(mu)
} # End of selection function
################# Main program #################
set.seed(100) # Initialize random number generator
N <- 100 # Set population size
MaxGen <- 2000 # Number of generations
Output <- matrix(0,MaxGen,2) # Create file for output
h2 <- 0.5 # Set heritability
Vp <- (.1)^2 # Set Phenotypic variance
Va <- Vp*h2 # Calculate Additive genetic variance
Ve <- Vp-Va # Calculate Environmental variance
mu <- 1.5 # Initial 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 environmental values using normal distri-
bution
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
# Calculate new mean genetic value by applying fitness criterion
mu <- SELECTION(X) # apply SELECTION
} # End of Igen loop
# Plot trajectory over generations
plot(Output[1:MaxGen,1], Output[1:MaxGen,2], type¼‘l’,
xlab¼‘Generation’, ylab¼‘Trait value’)
mean(Output[500:MaxGen,2]) # Mean phenotype
sd(Output[500:MaxGen,2]) # SD of mean
OUTPUT: (Figure 4.3)
GENETIC MODELS 247