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

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

[xx,yy] ¼ meshgrid(x, y)

zz ¼ x.^2 þ y.^2

The data are now in a format suitable for plotting the contours using contour or a

3-D plot using the R function persp (Figure 1.8) or surfc in MATLAB.

R CODE:

par(mfrow¼c(1,2)) # Divide graphics page into two

contour(x,y, z.matrix, xlab¼“x”, ylab¼“y”) # Contour plot

persp(x, y, z.matrix, theta¼30, phi¼10) # 3-D plot

The parameters theta and phi are angles deﬁning the viewing direction: theta

gives the azimuthal direction and phi the colatitude.

MATLAB CODE:

subplot(1,2,1); % Divide graphics page in two and plot contour in left

[C,h] ¼ contour(x, y, zz) % Create contour plot

02468

z.matrix

13 12

–10

–5 0 5 10

Figure 1.8 An example of contour and 3‐D plots.

48 MODELING EVOLUTION

% clabel(C,h) rotates the labels and inserts them in the contour

lines

clabel(C,h);

xlabel(‘x’); ylabel(‘y’); % Add text

subplot(1,2,2); % Divide graphics page in two and plot contour on

right

surfc(xx,yy,zz); % Plot a 3-D surface

xlabel(‘Foraging’); ylabel(‘Vigilance’); zlabel(‘Fitness’) % Add text

Now consider the analysis of variation in P.mig and P.surv in model 5. In this

case I want to plot variation in mean population size for the ranges

P.mig <- seq(from¼0.1, to¼0.9, length¼10) # Proportion migrating

P.surv <- seq(from¼0.8, to¼0.9, length¼10) # Survival rate of mi-

grants

To do this I make the previous main program into a function

MAIN.PROG <- function(D)

{

P.mig <- D[1]

P.surv <- D[2]

Same lines as previously except for deletion of plotting codes

} # End of function MAIN.PROG

Ihavelefttheset.seed code in MAIN.PROG which means that the same

sequence of ra ndom numbers is used for eac h comb ination: thus I replicate

exactly the same envir onmental variation for each simulation. Now I can use

expand.grid to generate the necessary variation and pass this to MAIN.

PROG. However, if I do this exactly as before I will get an error message,

because I need to pass only one combination at a time. To pass only one

combination per run use the R function apply telling the function to use

rows:

d <- expand.grid(P.mig,P.surv)

z <- apply(d, 1, MAIN.PROG)

z.matrix <- matrix(z, length(P.mig), length(P.surv))

par(mfrow¼c(1,2))

contour(P.mig, P.surv,z.matrix, xlab ¼“P.mig”, ylab¼“P.surv”)

persp(P.mig, P.surv, z.matrix, theta¼20, phi¼20)

OUTPUT: (Figure 1.9)

OVERVIEW 49

The results plotted in Figure 1.9 show that population size increases with

P.surv (not surprisingly) and that it is largest at an intermediate value of P.mig.

Step 20: Adding inheritance: functions pnorm, dnorm and numerical

integration

The observation from the above analysis of the potential importance of migration

in a heterogeneous environment raises the question of what is the optimal

migration rate. This question is made difﬁcult to answer analytically because

at the metapopulation level ﬁtness is frequency-dependent (Roff 1994a). Migrants

can be reasonably assumed to suffer an increased mortality rate (in models 4 and

5) and a decrease in reproductive opportunities by a lack of time and by

the allocation of energy to migration, both in terms of the energy used in migra-

tion and the energy sequested for the capability of migration, as found in wing-

dimorphic insects (Roff 1996; Roff and Fairbairn 2007). For the present analysis

I shall use a threshold model for the genetic basis of migrant and nonmigrant

phenotypes.

P.mi

P.surv

0.2

0.80 0.82 0.84 0.86

0.88 0.90

0.4 0.6 0.8

250

100

P.mig

P.surv

z.matrix

200

150

Figure 1.9 Contour and perspective plots for model 5.

50 MODELING EVOLUTION

1.4.7 Mathematical assumptions of model 6

1. Assumptions of model 5 hold for model 6 except that migration propensity is

inherited.

2. Mating is at random and occurs prior to migration.

3. Individuals can be divided into two classes, nonmigrants and potential mi-

grants. Individuals in the latter class migrate with a probability P

mig

and survive

the migratory episode with probability P

surv

. Additionally, potential migrants

suffer a reduction in reproductive ﬁtness because of their allocation of re-

sources to the capability of migration (e.g., presence or absence of a functional

ﬂight apparatus in wing-dimorphic insects).

4. Migratory propensity is inherited as a quantitative trait as speciﬁed by the

threshold model. According to this model there is a continuously normally

distributed underlying trait called the liability: individuals above a threshold

develop into one type of morph, while individuals below the threshold develop

into the alternate (Figure 1.10). Without loss of generality, we can set the

threshold at zero and the variance of the distribution at 1. Assuming, again

without loss of generality, that nonmigrants are those individuals lying above

the threshold the proportion of nonmigrants in the population is given by the

normal distribution function

P ¼



x  mÞdx; where fðxÞ¼

ﬃﬃﬃﬃﬃﬃ



ð1:36Þ

where m is the mean value of the liability. This function cannot be integrated

symbolically but there is an R function pnorm and a MATLAB function normcdf

that provide a solution. The integral from zero to inﬁnity, given a mean of Mu

(which can also be vector of values) is obtained as

P <- pnorm( 0, mean¼ -Mu, sd¼1) # R code

P ¼ normcdf(0, Mu, 1);% MATLAB code

Phenotype BPhenotype A

Liability

Threshold

Frequency

Figure 1.10 Graphical illustration of the threshold model. Individuals with liabilities below

the threshold develop into phenotype A, whereas individuals with liabilities above the

threshold develop into phenotype B.

OVERVIEW 51

The function (x) can be calculated using the R function dnorm or the MATLAB

function normcdf. Numerical integration can be done in R using the function

integrate and in MATLAB the functions int or quad (see, e.g., Scenario 4 of

chapter 2). The mean values of nonmigrants, X

, and potential migrants, X

can be calculated from the respective truncated normal distributions, giving

¼ m þ

fðmÞ

; X

¼ m 

fðmÞ

1  P

ð1:37Þ

Given random mating prior to migration, the phenotypic value of the liability of

nonmigrants, Y

, and potential migrants, Y

, is (Roff 1994)

¼ m þ

fðmÞh

; Y

¼ m 

fðmÞh

2ð1  PÞ

ð1:38Þ

which can be coded in R as (respectively)

Y.nonmigrants <-Muþ dnorm(0, mean¼Mu, sd ¼1)*h2/(2*P)

Y.migrants <- Mu - dnorm(0, mean¼Mu, sd¼1)*h2/(2*(1-P))

where h

is the heritability of liability (set at 0.5 in the present model). The mean

phenotypic value of offspring in the ith patch at generation tþ1, m

tþ1,i

is given by

(Roff 1994)

tþ1;i

t;i

½Y

NM;t;i

t;i

þ Y

PM;t;i

ð1  P

t;i

Þð1  P

mig

ÞCþ



t;i

½ð1  P

t;i

Þþð1  P

mig

ÞCþ

t;T

ð1:39Þ

where n is the number of patches, C is the “cost” incurred by potential migrants for

the ability to migrate (whether done so or not) and set to 0.6 in the present simula-

tion (this cost is consistent with the loss in fecundity of the winged morph of the

sand cricket [Roff 1984a]), N

t,i

is the number in the ith patch before migration, N

i,T

the total number of migrants, and Y



is the weighted phenotypic value of the

offspring from migrants. The latter two variables are deﬁned as

t;T

¼ CP

surv

mig

i¼1

t;i

ð1  P

t;i



¼ CP

surv

mig

i¼1

t;i

PM;t;i

ð1  P

t;i

ð1:40Þ

Equations (1.39) and (1.40) can be coded in R as

Emigrants <- P.mig*Npop*(1-P) # Nos of emigrants

Nos.migrants <- P.surv*sum(Emigrants) # N

t,T

Y.star <- P.mig*sum(Npop*P.surv*Y.migrants*(1-P))

Mu <- (Npop*(Y.nonmigrants*PþY.migrants*(1-P)*(1-P.

mig)*Cost) þ (Y.star/N.patches))/(Npop*(Pþ (1-

P)*(1-P.mig)*Cost) þ Nos.migrants/N.patches)

Finally, we compute the change in population size within each patch and return

the new population size and new mean liability for each population:

52 MODELING EVOLUTION

# New population size before reproduction

Npop <- Npop - Emigrants þ Nos.migrants/N.patches

Npop <- Npop*LAMBDA # Population size before

constraints

Npop[Npop<1] <- 0 # Check for extinction

Npop[Npop>1000] <- 1000 # Carrying capacity

# Return the vector of new population sizes and means

return (c(Npop,Mu))

The full coding is

R CODE:

rm(list¼ls()) #Clearmemory

# Population and inheritance function

POP <- function( MAX.Lambda,Npop, N.patches,Mu, P.surv, P.mig)

{

h2 <- 0.5; Cost <- 0.6 # parameters

LAMBDA <- runif(N.patches,min¼0, max¼MAX.Lambda) # randomlamb-

das

P <- pnorm( 0, mean¼ -Mu, sd¼1) # Proportion of nonmigrants

Y.nonmigrants <-Muþ dnorm(0,mean¼Mu, sd¼1)*h2/(2*P)

Y.migrants <- Mu - dnorm(0, mean¼Mu, sd¼1)*h2/(2*(1-P))

Emigrants <- P.mig*Npop*(1-P) # vector of surviving

emigrants

Nos.migrants <- P.surv*sum(Emigrants)# ¼N(t,T)

Y.star <- P.mig*sum(Npop*P.surv*Y.migrants*(1-P))

Mu <- (Npop*(Y.nonmigrants*PþY.migrants*(1-P)*(1-P.mig)

*Cost) þ (Y.star/N.patches))/(Npop*(Pþ (1-P)*(1-P.

mig)*Cost) þ Nos.migrants/N.patches)# Calculate Nos in

new populations

Npop <- Npop - Emigrants þ Nos.migrants/N.patches

Npop <

- Npop*LAMBDA # Population size before constraints

Npop[Npop<1]

0 # Check for extinction

Npop[Npop>1000] <- 1000 # Carrying capacity

return (c(Npop,Mu))# Return the vector of new popnsizesand means

} # End of function

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

set.seed(100) # set seed

MAXGEN <- 2000 # Set maximum number of generations

N.init <- 20 # Initial population size

MAX.LAMBDA <- 2.2 # Maximum value of lambda

N.patches <- 10 # Number of patches

P.surv <- 0.95 # Survival rate of migrants

P.mig <-.8 # Proportion of potential migrants migrating

Mu <- matrix(0,N.patches,1) # Initial mean liability values

Npop <- matrix(N.init, N.patches, 1) # Initialise populations

OVERVIEW 53

Npop.Sizes<- matrix(0,MAXGEN) # Pre-assign storage for means

Npop.Sizes[1] <- mean(Npop) # Store 1st generation

mean population size

N.extinct <- matrix(0,N.patches,1)# Assign storage for nos

extinct

# Pre-assign space for mean propn non-migrants

Mean.nonmig <- matrix(0,MAXGEN) # Storage for propn nonmigrants

# Mean propn nonmigrants

Mean.nonmig[1] <- mean(pnorm( 0, mean¼ -Mu, sd¼1))

Mean.mig <- matrix(0, MAXGEN) # Storage for propn migrants

Mean.mig[1] <- 1-Mean.nonmig[1] # Proportion migrants

Igen <- 1 # Initial generation number

while ( Igen<MAXGEN && Npop.Sizes[Igen]>0) # Enter while loop

{

Igen <- Igenþ1 # Increment generation counter

# Get new population sizes and mean liabilities

OUT <- POP(MAX.LAMBDA, Npop, N.patches, Mu, P.surv, P.mig)

Npop <- OUT[1:N.patches] # Vector of Population sizes

n1 <- N.patchesþ1; n2 <-2*N.patches # Range for mean liabilities

Mu <- OUT[n1:n2] # Mean liabilities

P.nonmigrants <- pnorm( 0, mean¼-Mu, sd¼1) # Vector of prop

nonmigrants

# Mean proportion of nonmigrants in metapopulation

Mean.nonmig[Igen] <- sum(Npop*P.nonmigrants)/sum(Npop)

# mean proportion of population migrating

Mean.mig[Igen] <- sum(Npop*(1-P.nonmigrants)*P.mig)/sum(Npop)

Npop.Sizes[Igen] <- mean(Npop) # Store mean population size

N.extinct[Igen] <- length( Npop [Npop ¼¼0]) # Store nos of ex-

tinct popns

}

par(mfcol¼c(2,2)) # Divide graphics page into four quadrants

Gen <- seq(1,Igen) # vector of generation numbers

plot(Gen, Npop.Sizes[1:Igen], xlab¼“Generation”, ylab=“Mean

populationsize”, type¼“l”) # Mean population size over generations

plot(Gen, N.extinct[1:Igen], xlab¼“Generation”, ylab

¼“Number

pops extinct”, type¼“l”, ylim¼c(0,N.patches)) # Nos extinct

over generations

plot(Gen, Mean.nonmig[1:Igen], xlab¼“Generation”, ylab¼“Mean

Proportion of nonmigrants”, type¼‘l’) # Mean proportion of nonnon-

migrants over generations

plot(Gen, Mean.mig[1:Igen], xlab¼“Generation”, ylab¼“Mean

Proportion of migrants”, type¼‘l’) # Mean proportion of migrants

over generations

OUTPUT: (Figure 1.11)

54 MODELING EVOLUTION

In addition to the population size data the program also plots the mean propor-

tion of nonmigrants and the mean proportion of the population that actually

migrate. As expected from the contour analysis, the system evolves to an interme-

diate level of migration. An interesting question, which I leave the reader to

address is whether or not the evolutionarily stable proportion maximizes popula-

tion size.

1.5 Summary of modeling approaches descr ibed in this book

1.5.1 Fisherian optimality analysis (Chapter 2)

Fisher’s general analysis of evolution was based on the characteristic equation and

the maximization of the Malthusian parameter, r (see Section 2.2). This approach

has formed the backbone of much of the study of the evolution of trait variation.

While I do not wish to imply that Fisher used only this approach, I think that he

can be acknowledged as its originator. To distinguish this type of analysis from

others in the book, which are also optimality models in one sense or another,

I shall refer to it as “Fisherian” optimality analysis. The general assumptions of

these models are

Generation

0246810

0.50

0.15 0.25 0.35 0.45

0.60 0.70 0.80

200 400 600 800

500 1,000 1,500 2,000 0 500 1,000 1,500 2,000

0 500 1,000 1,500 2,000

Mean Proportion of migrants Mean Proportion of nonmigrants

Number of pops extinct

Mean population size

Figure 1.11 Mean population size, mean proportion of nonmigrants and actual migrants,

and number of extinct populations for model 6. For display purposes the graph page has

been split into four quadrats using par(mfcol=c(2,2)).

OVERVIEW 55

1. Stable-age distribution.

2. If density-dependence occurs it does not affect the trait under study.

3. Selection is frequency-independent.

4. In a stable environment ﬁtness is maximized by r, R

or a ﬁtness component

whose maximization also maximizes one of the aforementioned ﬁtness mea-

sures. The way these ﬁtness measures are calculated depends upon the type of

environment.

5. In a temporally variable environment the appropriate operational measure of

ﬁtness is the geometric mean of r. With age structure the appropriate ﬁtness

measure is the dominant Lyapunov exponent and is most readily analyzed

using invasibility analysis (see Chapter 3).

6. Spatial variation can be handled by this approach using the framework devel-

oped simultaneously by Houston and McNamara (1992) and Kawecki and

Stearns (1993), as described by equations (1.28) and (1.29). This assumes an

equilibrium population.

Examples: Age at ﬁrst reproduction and allocation to reproduction.

1.5.2 Invasibility analysis (Chapter 3)

Invasibility analysis is a technique that determines if a clone showing an alternate

suite of characteristics (generally differences in parameter values) can invade a

population. While it can be used for unstructured populations its strength lies

primarily in addressing the evolution of traits in structured populations in which

density-dependence occurs. The general properties are

1. The population is age-structured or stage-structured. Typically, the model is

written in matrix form (e.g., a Leslie matrix).

2. The population may be at equilibrium or show variable behavior (e.g., cyclic or

chaotic). If equilibrium can be assumed, then the model may be better analyzed

using the methods for “Fisherian” optimality models.

3. Density-dependence may occur and may be an important selective factor.

Examples: Evolution of parameters in density-dependent functions (e.g., the

Ricker function) and delayed maturity in a stage-structured model.

1.5.3 Genetic models (Chapter 4)

The preceding methods of analysis presume that the genetic mode of inheritance

is not an impediment to the population evolving to the optimal combination of

trait values. “Fisherian” optimality analysis focuses upon the optimal combination

and has little to say about the trajectory or time course from one combination to

another. Invasibility analysis can follow the change in frequency of mutant clones

but whether this is relevant to other modes of inheritance is not addressed (in fact,

as shown in Chapter 3, the genetic mode of inheritance can profoundly inﬂuence

56 MODELING EVOLUTION

the outcome of evolution). In Chapter 3 I consider three types of genetic models,

one focused at the population level and the other two focused on individual-based

models. Potential advantages of genetic models are

1. The evolutionary trajectory can be followed.

2. Variation about the equilibrium combination can be assessed.

3. Density-dependence and frequency-dependence are easily incorporated into

the individual-based approaches (examples of the latter are given in Chapter 6).

4. Age structure is readily accommodated (though bookkeeping can be tedious).

5. The model is itself a simple numerical method of locating the optimal combi-

nation. Used in this sense the actual genetic parameters may be irrelevant and

can be selected to minimize the “search” time (e.g., set genetic correlations to

zero). In rare cases parameter values may alter the equilibrium combination,

but this can be easily investigated by varying parameter values.

Examples: Evolution of multiple traits in density-dependent populations and trait

variation under the impact of stabilizing or directional selection.

1.5.4 Game theoretic models (Chapter 5)

Game theoretic models deal with interactions among individuals, typically be-

tween pairs of individuals. Such interactions, or “games,” may be frequency-

independent or frequency-dependent. In both cases ﬁtness is deﬁned by a payoff

matrix that gives the increment (or decrement) in ﬁtness associated with all

possible pairs of interactions. The frequency-independent case is really a subset

of “Fisherian” optimality analysis and most game theoretic models involve fre-

quency-dependent games. General characteristics of game theoretic models are

1. Typically, population dynamics is not an important factor, a stable population

being a frequent, but unspeciﬁed, assumption.

2. In simple models analytical solutions are possible but in more complex cases an

individual-based models are used.

3. Clonal or explicit genetic models can be used.

Examples: Hawk-Dove game and territorial–satellite behavior.

1.5.5 Dynamic programming (Chapter 6)

Individuals are continually faced with decisions that affect their ﬁtness. For

example, a parasitic wasp must make a decision (I do not imply that this is a

conscious act) on how many eggs to lay in a host: Too many larvae coul d compete

with each other and hence have a reduction in ﬁtness, but the female might not

ﬁnd another host and so laying too few larvae in the present host could result in a

reduced ﬁtness. In some cases the optimal decision set can be estimated using

one of the previous approaches but in many cases the sequence of decisions to be

made during the life of an organism is so complex that these approaches are

OVERVIEW 57