Roff D.A. Modeling Evolution: An Introduction to Numerical Methods

Подождите немного. Документ загружается.

The set of combinations of R and S are generated using the R function expand.

grid. Because we require the storage of three different values it is most conve-

nient to iterate over these combinations rather than using the apply function.

Combinations that are not physically possible (total proportion exceeding 1) are

skipped. The three resulting vectors are converted into a matrix for contour

plotting and the three zero isoclines of allele frequencies (use levels=0 in call

to contour function) are plotted on the same plot (after the ﬁrst plot add draw-

labels=FALSE to the call to the contour function). The calculation of the con-

tours for the proportions of rocks, papers, and scissors is a little trickier because

we have to convert the allele frequencies into morph frequencies. To do this we

make use of the R function contourLines, which calculates the set of x,y

coordinates for the user-speciﬁed contour lines. In this case we want only the

zero contours (so pass levels=0). First we obtain the relevant coordinates for the

matrix X (say), which has dimensions Proportion by Proportion:

Data <contourLines(Proportion, Proportion, X, levels=0)

The object Data is a list from which we must extract the relevant vectors. This is

done using the R function unlist. An added complication is that contourLines

might construct several zero isocline combinations: each set has to be extracted

separately. Once this is done the isocline can be plotted using the x,y coordinates.

This is done in the user-supplied function ISOCLINE.

R CODE:

rm(list=ls()) # Remove all objects from memory

# Function to calculate new proportions

FITNESS <- function(Morph, PayoffMatrix, Genotypes)

{

# To take account of integer values Npop is recalculated here

Npop <- sum(Genotypes)

# Match individuals up to find fitness for each

# Create a randomized vector of opponents

Opponent <- sample(Morph)

Fitness <- matrix(0,Npop,1) # Allocate space for Fitness vector

# Iterate over the Payoff matrix

for (Receiver in 1:3 ) # Ind receiving payoff 1=Rock 2 ¼ Scissors

3=Paper

{

for (I.Opponent in 1:3) # Opponent 1=Rock 2 ¼ Scissors 3=Paper

{

Fitness[Morph==Receiver & Opponent==

I.Opponent]<-5þ PayoffMatrix[Receiver,I.Opponent]

}}

# Calculate ranges for genotypes to count alleles

n0 <- Genotypes[1] # Nos of RR. This is for completness

n1 <-n0þ1 # Starting row of RS

308 MODELING EVOLUTION

n2 <-n1þGenotypes[2]-1 # Ending row of RS

n3 <-n2þ1 # Starting row of RP

n4 <-n3þGenotypes[3]-1 # Ending row of RP

# Sum R alleles RR RS RP

R.alleles <- 2*sum(Fitness[1:n0])þsum(Fitness[n1:n2])þsum

(Fitness[n3:n4])

# Number of S alleles

R.alleles <- 2*sum(Fitness[1:n0])þsum(Fitness[n1:n2])þsum

(Fitness[n3:n4])

# Number of S alleles

n5 <-n4þ1 # Starting row of SS

n6 <-n5þGenotypes[4]-1 # Ending row of SS

n7 <-n6þ1 # Staring row of SP

n8 <-n7þGenotypes[5]-1 # Ending row of SP

# Sum S alleles SS RS SP

S.alleles <- 2*sum(Fitness[n5:n6])þsum(Fitness[n1:n2])þsum

(Fitness[n7:n8])

# Number of P alleles

n9 <-n8þ1 # Starting row of PP

# Number of P alleles PP RP SP

P.alleles <- 2*sum(Fitness[n9:Npop])þsum(Fitness[n3:n4])þ

sum(Fitness[n7:n8])

# Proportion of each allele

Prop.R <- R.alleles/(R.alleles þ S.alleles þ P.alleles) # Propn R

allele

Prop.S <- S.alleles /(R.a lleles þ S.alleles þ P.alleles) # Propn S

allele

Prop.P <- P.alleles/(R.alleles þ S.alleles þ P.alleles) # Propn P

allele

turn(c(Prop.R, Prop.S, Prop.P)) # Returnproportion of alleles

} # End of function

###################### FUNCTION GENOTYPE ######################

GENOTYPE <- function(Prop.R, Prop.S, Prop.P, Npop)

{

# Calculate the genotypes

# First all genotypes that have Rock phenotype

RR <- round(Prop.R^2*Npop) # Number of RR genotypes

RS <- round(2*Prop.R*Prop.S*Npop) # Number of RS genotypes

RP <- round(2*Prop.R*Prop.P*Npop) # Number of RP genotypes

# Genotypes that have Scissors phenotype

SS <- round(Prop.S^2*Npop) # Number of SS genotypes

SP <- round(2*Prop.S*Prop.P*Npop) # Number of SP genotypes

# Genotypes that have Paper phenotype

PP <- round(Prop.P^2*Npop) # Number of PP genotypes

N <-RRþRSþRP þSSþSPþ PP # To account for rounding effects

GAME THEORETIC MODELS 309

P.Rocks <- (RRþRSþRP)/N # Proportion of Rocks

P.Scissors <- (SSþSP)/N # Proportion of Scissors

P.Paper <- PP/N # Proportion of Paper

return( c(RR, RS, RP, SS, SP, PP, P.Rocks, P.Scissors, P.Paper))

} # End of function

###################### FUNCTION ISOCLINE ######################

ISOCLINE <- function(X, ADD, LTY) # Function to plot zero isoclines

{

# X is the matrix of changes in proportion

# ADD is a flag. ADD=0 tells function to start a new plot

# LTY is the type of line to be drawn

# Call R function contourLines to get xy coordinates of zero isoclines

Data <- contourLines(Proportion, Proportion, X, levels=0)

b <- data.frame(unlist(Data)) # Unlist Data and convert to a

data frame

N.zeros <- length(b[b==0]) # Number of zero isoclines

for ( i in 1:N.zeros) # Iterate over isoclines

{

Prop.R <- Data[[i]]$x # Proportion R allele

Prop.S <- Data[[i]]$y # Proportion S allele

Prop.P <- 1- Prop.R - Prop.S # Proportion P allele

RR <- Prop.R^2 # Frequency of RR genotypes

RS <- 2*Prop.R*Prop.S # Frequency of RS genotypes

RP <- 2*Prop.R*Prop.P # Frequencyof RP genotypes

P.Rocks <-RRþRSþRP # Proportion of Rocks

SS <- Prop.S^ 2 # Frequency of SS genotypes

SP <- 2*Prop.S*Prop.P # Frequency of SP genotypes

P.Scissors <-SSþSP # Proportion of Scissors

# Plot lines (Have to check if new plot is requested)

if(i==1 && ADD==0) {plot(P.Rocks,P.Scissors, type=’l’,

xlab=“Proportion of Rock”, ylab=“Proportion of Scissors”,

xlim=c(0,1), ylim=c(0,1),lty=LTY)}

else{ lines(P.Rocks,P.Scissors, lty=LTY) }

}

} # end of function

##################### MAIN PROGRAM #####################

set.seed(100) # Initialize random number generator

Npop <- 2000 # Set population size

# Create a sequence of proportions for the two R and S alleles

Nos.P <- 30 # Number of divisions

Proportion <- seq(from=0.01, to=0.95, length=Nos.P)

310 MODELING EVOLUTION

# Use expand.grid to generate all possible combinations

RxS <- expand.grid(Proportion, Proportion)

Data.Rocks <- matrix(9,Nos.P*Nos.P,1) # Allocate space for R

output

Data.Scissors <- Data.Rocks # Allocate space for S output

Data.Paper <- Data.Rocks # Allocate space for P output

Morph <- matrix(0,Npop,1) # Allocate space for Phenotype

vector

# Set up fitness matrix. Note column-wise fill

Epsilon <- 0.1 # Deficit when the same morph types interact

PayoffMatrix <- matrix(c(-Epsilon,-1,1,1,-Epsilon,-1,-1,1,

-Epsilon),3,3)

# Iterate over all possible combinations

Total.Combinations <- Nos.P*Nos.P # Total numberof combinations

for (Ith.comb in 1:Total.Combinations)

{

Prop.R <- RxS[Ith.comb,1] # Proportion of R allele

Prop.S <- RxS[Ith.comb,2] # Proportion of S allele

Prop.P <- 1-Prop.R-Prop.S # Proportion of P allele

Total <- Prop.RþProp.SþProp.P # Sum of proportions

if(Total > 0) # Check that combination is permissable

{

D1 <- GENOTYPE(Prop.R, Prop.S, Prop.P, Npop) # Calculate the

genotypes

Genotypes <- D1[1:6] # Numbers of each genotype

Prop.Morphs <- D1[7:9] # Proportions of each morph

Morph[1:Npop] <- 3 # Set initially to Paper

Nos.of.Rocks <- sum(Genotypes[1:3]) # RRþRSþRP

Nos.of.Scissors <- sum(Genotypes[4:5]) # SSþSP

Nos.of.Papers <- Genotypes[6] # PP genotype

Morph[1:Nos.of.Rocks] <- 1 # Set rows 1 to Nos.of.Rocks to Rocks

n1 <- Nos.of.Rocksþ1 # Starting row for Scissors

n2 <-n1þNos.of.Scissors-1 # Ending row for Scissors

Morph[n1:n2] <- 2 # Set these rows to Scissors (¼2)

# Calculate new proportions by applying fitness criterion

Propns <

- FITNESS(Morph,PayoffMatrix,Genotypes)

D2 <-

GENOTYPE(Propns[1], Propns[2], Propns[3], Npop)

Data.Rocks[Ith.comb] <- D1[7]-D2[7] # Change in proportion of

Rocks

Data.Scissors[Ith.comb] <- D1[8]-D2[8] # Change in propor-

tion of Scissors

Data.Paper[Ith.comb] <- D1[9]-D2[9] # Change in propor-

tion of Paper

GAME THEORETIC MODELS 311

} # End of if statement

} # End of Ith.comb loop

# Convert vectors to matrices Note that fill is by column not row

Delta.Rocks <- matrix(Data.Rocks, Nos.P, Nos.P, byrow=F)

Delta.Scissors <- matrix(Data.Scissors, Nos.P, Nos.P, byrow=F)

Delta.Paper <- matrix(Data.Paper, Nos.P, Nos.P, byrow=F)

par(mfrow=c(2,2)) # 4 plots per page

# Plot zero isoclines for alellic frequencies

contour(Proportion, Proportion, Delta.Rocks, xlab=“Proportion

of R allele”, ylab=“Proportion of S allele”, levels=0, drawla-

bels=FALSE)

contour(Proportion, Proportion, Delta.Scissors,lty=2,

levels=0, add=TRUE, drawlabels=FALSE)

contour(Proportion, Proportion, Delta.Paper, lty=3, levels=0,

add=TRUE, drawlabels=FALSE)

# Get zero isoclines for proportion of each morph

ISOCLINE(Delta.Rocks, 0, 1) # Plot Rock zero isocline

ISOCLINE(Delta.Scissors, 1, 2) # Plot Scissors zero isocline

ISOCLINE(Delta.Paper, 1, 3) # Plot Paper zero isocline

points(0.333,0.333, cex=2) # Add predicted ESS. cex ¼ large

symbol

OUTPUT: (Figure 5.8)

The graphical output shows that an equilibrium is possible. Note that the allele

frequencies differ from each other, as expected from the genotype–phenotype

map.

Proportion of R allele

0.0 0.2 0.4 0.6 0.8

Proportion of Rock

Proportion of Scissors

0.0

0.0 0.2 0.4 0.6 0.8 1.0

Proportion of S allele

0.0 0.2 0.4 0.6 0.8

0.2 0.4 0.6 0.8 1.0

Figure 5.8 Zero isoclines for the R‐P‐S Mendelian model (Scenario 6): solid line, rock;

dashed line, scissors; and dotted line, paper. The broken line at 45° can be ignored.

The predicted ESS is indicated by the circle in the right‐hand plot.

312 MODELING EVOLUTION

5.8.4 Finding the ESS using a numerical approach

R CODE:

rm(list=ls()) # Remove all objects from memory

# Function to calculate new proportions

FITNESS <- function(Morph,PayoffMatrix,Npop,Genotypes)

{

SAME AS PREVIOUS CODING

}

GENOTYPE <- function(Prop.R, Prop.S, Prop.P, Npop)

{

SAME AS PREVIOUS CODING

}

##################### Main program #####################

set.seed(100) # Initialize random number generator

Npop <- 2000 # Set population size

MaxGen <- 2000 ) # Number of generations

Output <- matrix(0,MaxGen,7) # Create file for output

Prop.R <- 0.33 # Proportion Rocks

Prop.S <- 0.33 # Proportion Scissors

Prop.P <- 1-Prop.R-Prop.S # Proportion Papers

Morph <- matrix(0,Npop,1) # Allocate space for Phenotype

vector

# Set up fitness matrix. Note column-wise fill

Epsilon <- 0.1 # Deficit when the same morph types interact

PayoffMatrix <- matrix(c(-Epsilon,-1,1,1,-Epsilon,-1,-1,1,

-Epsilon),3,3)

for (Igen in 1:MaxGen) # Iterate over generations

{

D1 <- GENOTYPE(Prop.R, Prop.S, Prop.P, Npop) # Calculate the

genotypes

Genotypes <- D1[1:6] # Numbers of each genotype

Prop.Morphs <- D1[7:9] # Proportions of each morph

Morph[1:Npop] <- 3 # Set initially to Paper

Nos.of.Rocks <- sum(Genotypes[1:3]) # RRþRSþRP

Nos.of.Scissors <- sum(Genotypes[4:5]) # SSþSP

Nos.of.Papers <- Genotypes[6] # PP genotype

Morph[1:Nos.of.Rocks] <- 1 # Set rows 1 to Nos.of.Rocks to Rocks

n1 <- Nos.of.Rocksþ1 # Starting row for Scissors

n2 <-n1þNos.of.Scissors-1 # Ending row for Scissors

Morph[n1:n2] <- 2 # Set these rows to Scissors(¼2)

Output[Igen,1] <- Igen # Store Generation number

N <- sum(Genotypes) # To avoid rounding problems

GAME THEORETIC MODELS 313

Output[Igen,2] <- Nos.of.Rocks/N # Proportion of Rocks

Output[Igen,3] <- Nos.of.Scissors/N # Proportion of Scissors

Output[Igen,4] <- Nos.of.Papers/N # Proportion of Papers

Output[Igen,5] <- Prop.R # Frequency of allele R

Output[Igen,6] <- Prop.S # Frequency of allele S

Output[Igen,7] <- Prop.P # Frequency of allele P

# Calculate new proportion of Rocks by applying fitness criterion

Propns <- FITNESS(Morph,PayoffMatrix,Genotypes)

Prop.R <- Propns[1] # Frequency of allele R

Prop.S <- Propns[2] # Frequency of allele R

Prop.P <- Propns[3] # Frequency of allele R

} # End of Igen loop

par(mfrow=c(2,2)) # 4 plots per page

plot(Output[,1], Output[,2], type=’l’, xlab=’Generation’,

ylab=’Morph Proportions þ X’, ylim=c(0.0,1.8)) # Plot Rocks

lines( Output[,1], Output[,3]þ0.5, lty=2) # Add Scissors to plot

lines( Output[,1], Output[,4]þ1.0, lty=3) # Add Papers to plot

# Add predicted line to plots

lines(Output[,1], rep(0.3333,MaxGen)) # “Predicted” ESS

lines(Output[,1], rep(0.3333þ0.5,MaxGen)) # “Predicted” ESS

lines(Output[,1], rep(0.3333þ1,MaxGen)) # “Predicted” ESS

# Plot Allele frquencies

plot(Output[,1], Output[,5], type=’l’, xlab=’Generation’,

ylab=’Allele Proportions’, ylim=c(0.0,1)) # Plot Rocks

lines( Output[,1], Output[,6], lty=2) # Add Scissors to plot

lines( Output[,1], Output[,7], lty=3) # Add Papers to plot

# Phase plots showing proportion of two morphs

plot(Output[,2],Output[,3],type=’l’, xlab=’Rocks’, ylab=’S-

cissors’)

plot(Output[,2],Output[,4],type=’l’, xlab=’Rocks’, ylab=’Pa-

pers’)

# plot(Output[,3],Output[,4],type=’l’, xlab=’Scissors’, ylab=’

Papers’)

# Print mean proportions starting at generation 400

print(’ Mean proportions (R,P,S) from Generation 400 to MaxGen’)

c(mean(Output[400:MaxGen,2]),mean(Output[400:MaxGen,3]),

mean(Output[400:MaxGen,4]))

print(’ Mean allele freqs (R,P,S) from Generation 400 to MaxGen’)

c(mean(Output[400:MaxGen,5]),mean(Output[400:MaxGen,6]),

mean(Output[400:MaxGen,7]))

OUTPUT: (Figure 5.9)

314 MODELING EVOLUTION

[1] “Mean proportions (R,P,S) from Generation 400 to MaxGen”

[1] 0.33980440.32908460.3311111

[1] “Mean allele freqs (R,P,S) from Generation 400 to MaxGen”

[1] 0.18938080.23929430.5713249

Both the allelic and morph frequencies show damped cyclic trajectories, with

means being as indicated from the earlier graphical analysis (see Figure 5.9). The

proportion of each morph oscillates about 1/3 but the allelic frequencies ﬂuctuate

about 0.19 (R), 0.24 (S), and 0.57 (P).

5.9 Scenario 7: Rock-Paper-Scissors: a quantitative genetic

model

Although one might imagine this to be a more difﬁcult model to code than

the Mendelian model it actually turns out to be somewhat simpler. As

with the Mendelian model, there are a number of ways in which genotype can

Rocks

0.1

0.1 0.2 0.3

Scissors

0.4 0.5 0.6

0.1 0.2 0.3

Papers

0.4 0.5 0.6

0.2 0.3 0.4 0.5 0.6

Rocks

Generation

0.1

0 500 1,000 1,500 2,000

Generation

0.0

0.0 0.2 0.4 0.6 0.8 1.0

Morph Proportions + X

Allele Proportions

0.5 1.0 1.5

500 1,000 1,500 2,000

0.2 0.3 0.4 0.5 0.6

Figure 5.9 Output from the Mendelian model of the R‐P‐S game. The traces for the

three proportions over time have been separated by adding 0.5 to the scissors

(middle trace) and 1.0 to papers (top trace). Predicted ESS values are shown as the

horizontal lines (ESS = 0.333 but increased for display purposes by 0.5 and 1.0 for

scissors and papers, respectively). Allele proportions have not been separated. The two

phase plots start at generation zero.

GAME THEORETIC MODELS 315

translate into phenotype. I shall assume here a simple model in which the three

morphs are determined by three different traits that are uncorrelated.

5.9.1 General assumptions

These are the same as that of Scenario 5.

5.9.2 Mathematical assumptions

1. The payoff matrix is as shown in Table 5.3. In this case when

¼ 0.1 the model

was completely unstable, one of the morphs eventually reaching 100%, which

one depending on chance (i.e., initial starting conditions such as the random

number generator). More interesting behavior was found with

¼ 0.5 and so

this is what is used here. However, it is important to note that analyses need to

consider a range of values in order to arrive at generalities.

2. Population size is ﬁnite and constant.

3. Only one interaction occurs per individual.

4. Morph is determined by a polygenetic system in the following manner:

There are three uncorrelated traits, each of which is a standard quantitative

trait (i.e., normally distributed). The trait R

underlies the rock morph, the trait

underlies the scissors morph, and the trait P

underlies the paper morph.

The morph that is phenotypically expressed is that with the highest value:

a plausible biological scenario for this could be the action of three hormones.

5. Mating is at random and genotypes at each generation are determined from

the multivariate normal distribution of genetic and environmental values (see

Chapter 5 for details of the quantitative genetic model). As with the previous

model, for simplicity, males and females are not distinguished, meaning that

the ﬁtnesses apply equally to each sex. In the lizard case it might be more

proper to assume that the R-P-S ﬁtnesses apply to the males and that the female

allele ﬁtnesses are equal. This would require separate coding for males and

females (see the Hawk-Dove model for how this can be implemented).

6. The mean genetic value of the next generation is determined from the weight-

ed (by ﬁtness) average of the previous genetic values. Thus the mean genetic

value of rock,



U;A

,is



ðt þ 1Þ¼

Pop

i¼1

ðtÞG

ðtÞ

Pop

i¼1

ðtÞ

ð5:21Þ

where W

(t) is the ﬁtness of the ith phenotype at generation t, G

(t) is the genetic

value of the ith phenotype at generation t, and N

Pop

is the population size.

5.9.3 A graphical analysis

A graphical analysis is not worthwhile in this case. There are three traits which

can vary, in principle, over an inﬁnite range. Equilibrium is expected whenever

316 MODELING EVOLUTION

the means of all three traits are the same. There is no simple way to graphically

display the result varying over the three parameter space. Therefore, we proceed

directly to the analysis over time.

5.9.4 Finding the ESS using a numerical approach

The initial variance–covariance matrices are determined by ﬁrst creating the

matrix H2 that contains the heritability of the three traits (here 0.5) along the

diagonal, the correlations (here 0) on the off-diagonals:

H2 ¼

0:50 0

00:50

000:5

ð5:22Þ

which in R code is given by

H2 <- matrix(c(0.5, 0, 0, 0, 0.5, 0, 0, 0, 0.5), 3,3)

Without loss of generality, the phenotypic variances are set at 1. The mean

environmental values are, by deﬁnition, equal to zero throughout the

simulation. Initially, the mean additive genetic values are set to zero; these

change as a result of selection. The diagonal elements of the covariance matrices

are the variances: hence for the phenotypic variance–covariance matrix, CovP,we

have

diag(CovP)<- c(vars[1], vars[2], vars[3])

where vars[1]< vars[2]< vars[3]< 1.

The additive genetic variances, V

, are obtained from the formula V

¼ V

where V

is the phenotypic variance and h

is the heritability. The

environmental variance, V

, is the difference between the phenotypic and

genetic variances V

¼ V

 V

. It is necessary to put a check in the coding since

must be positive. In R we can code this as

for ( i in 1:NX) # Iterate over traits

{

CovA[i,i] <- CovP[i,i]*H2[i,i] # ¼ Vp*h2

CovE[i,i] <- CovP[i,i]-CovA[i,i] # Environmental Variance

if(CovE[i,i] < 0) stop (print(c(“CovE cannot be”,i,j,CovE[i,i])))

} # end of i 1:NX loop

The covariances (off-diagonal elements) of the phenotypic and genetic covariance

matrices are obtained from the relationship Cov

¼ r

, where X and Y are

the two traits and r is the correlation. The environmental covariance is obtained

by subtraction: coding in R is

N.minus.1 <- NX-1

for( i in 1:N.minus.1)

{

jj <-iþ1

for(j in jj:NX)

GAME THEORETIC MODELS 317