Roff D.A. Modeling Evolution: An Introduction to Numerical Methods
Подождите немного. Документ загружается.
diag(Dummy) <- S[1:M] # Assign survivals to diagonal
Leslie.matrix <- matrix(0, Maxage, Maxage) # Pre-assign space
Leslie.matrix[1,] <- Fertility # Add fertilities to top row
Leslie.matrix[2:Maxage,1:M] <- Dummy # Add dummy to appropriate
space
# An alternate approach using a loop is shown below
#j <- 0; for (i in 2:Maxage){j <-jþ1 ;Leslie.matrix[i,j] <-S
[i-1]}
Eigen.data <- eigen(Leslie.matrix) # Call eigen
Lambda <- Eigen.data$values[1] # Get first eigenvalue
return(abs(log(Lambda))) # Return r
}
Optimum <- optimize(RCALC, interval¼ c(1,3), maximum¼TRUE)
Best.X <- Optimum$maximum # Optimum body size
Best.r <- Optimum$objective # Maximum r
# Print out results to 6 significant digits
print(c(’Best x ¼ ’, signif(Best.X, 6), ’Best r ¼ ’, signif
(Best.r,6)))
# Plot r.est vs x
n <- 50 # Nos of increments
x <- matrix(seq(from¼1, to¼2.5, length¼n)) # Values of x
r.est <- apply(x,1,RCALC) # Get values of r
plot(x, r.est,xlab¼“Body size, x”, ylab¼“r.est”, type¼’l’)
points( Best.X, Best.r, cex¼2) # Add point to graph at optimum
OUTPUT:
Figure not shown but same as Figure 2.5, except for added point at optimum.
[1] “Best x ¼ ” “1.93751” “Best r ¼ ” “1.49699”
which agrees, as expected, with the previous results. The matrix approach is
somewhat simpler in its coding compared to the summation approach of
the last section but not to that using the explicit function that relates r to x
(equation (3.25)).
3.4 Scenario 2: Adding density-dependence
We continue with the previous model but add density-dependence and use pair-
wise invasibility and elasticity analyses to locate the optimum body size.
3.4.1 General assumptions
1. All assumptions given in Scenario 1 hold.
2. Population size is limited by density-dependence.
188 MODELING EVOLUTION
3.4.2 Mathematical assumptions
1. All assumptions given in Scenario 1 hold.
2. Population size is controlled by a Ricker density-dependent function:
F
t
¼ F
t
ae
bN
t
ð3:26Þ
where N
t
is total population size, F
i
is the density-independent component of
fertility as defined by equation (3.20), and F
t
* is the density-dependent fertility.
3.4.3 Solving using R
0
as the fitness measure
Because the density-dependence does not directly affect the trait under consider-
ation the appropriate measure of fitness is not r but R
0
. The relevant equation is
equation (3.24) rewritten as
ða
F
þ b
F
xÞe
ða
S
þb
S
xÞ
1 e
ða
S
þb
S
xÞ
¼ R
0
ð3:27Þ
We can use the previous coding, modified for the change in function, to find the
value of x that maximizes R
0
(alternatively, one could find x such that
dR
0
dx
¼ 0).
R CODE:
rm(list¼ls()) # Clear workspace
RCALC <- function(x) # Function to calculate r
{
Af <-0; Bf<-16; As<-1; Bs<-0.5 # parameter values
A <-AfþBf*x; B <-AsþBs*x # For convenience
R0 <- A*exp(-B)/(1-exp(-B)) # R0
return(R0) # return value
}
# Call optimize to find best x
optimize(f¼RCALC,interval¼ c(0.1,3), maximum¼TRUE)$maximum
OUTPUT:
[1] 1.682795
Although the density-dependence does not directly involve body size it changes
the operational fitness measure and hence also the optimal body size. We now
examine the approaches of invasibility and elasticity analyses.
3.4.4 Pairwise invasibility analysis
The program follows the pattern outlined in the introduction with minor changes.
A general description of the functions follows:
1. LESLIE <- function(x,Maxage): This is the same function as in the previ-
ous scenario and constructs the Leslie matrix from the relevant equations.
INVASIBILITY ANALYSIS 189
2. DD.FUNCTION <- function(ALPHA, BETA, Fi, n) {Fi*ALPHA*exp
(-BETA*n)}: This passes the two density dependent parameters, the DI (densi-
ty-independent) fertility coefficient and total population size, and passes back
the new fertility as defined by equation (3.26). Parameter values are set at a ¼ 1, b
¼ 2 10
5
, which produces a stable equilibrium of 80782.
3. POP.DYNAMICS <- function(X): X contains the body size of the resident and
the body size of invader (this differs from the example in the introduction
which passed the multiplier for the invader). The function calculates the
growth rate of the invader. Unlike the example given in the introduction the
trajectory of the resident-only population is not followed.
4. Main program: This follows the approach outlined in the introduction. The
body size obtained from the elasticity analysis that follows is plotted onto the
contour surface.
R CODE:
rm(list¼ls()) # Clear workspace
LESLIE <- function(x,Maxage) # Function to generate Leslie matrix
{
M <- Maxage-1 # 1 less than the maximum age
age <- seq( from¼1, to¼Maxage) # vector of ages
Af <- 0 ;Bf <- 16; As <-1 ; Bs <- 0.5 # Parameter values
m <- rep(AfþBf*x, times¼Maxage) # number of female offspring
l <- exp(-(AsþBs*x)*age) # Survival to age
S <- matrix(0,Maxage,1) # Space for age-specific survival
S[1] <- l[1] # Survival to age 1
# Calculate the survival from t to tþ1
for ( i in 2:M) {S[i] <- l[i]/l[i-1]}
Fertility <- m*S # Top row of Leslie matrix
Dummy <- matrix(0,M,M) # Create a temporary matrix
diag(Dummy) <- S[1:M] # Assign survivals to diagonal
Leslie.matrix <- matrix(0, Maxage, Maxage) # Pre-assign space
Leslie.matrix[1,] <- Fertility # Add fertilities to top row
Leslie.matrix[2:Maxage,1:M] <- Dummy # Add dummy to appropriate
space
return(Leslie.matrix)
} # End of Leslie function
############### Density-dependence function ###############
DD.FUNCTION <- function(ALPHA,BETA,Fi,n) {Fi*ALPHA*exp
(-BETA*n)}
###############Population dynamics function ###############
POP.DYNAMICS <- function(X)
{
X.Resident <- X[1] # Body size of resident population
X.invader <
- X[2] # Body size of invader
ALPHA <-
1 ; BETA <- 2*10^-5 # Density dependence parameters
Maxage <- 50 # Maximum age
190 MODELING EVOLUTION
Resident.matrix <- LESLIE(X.Resident, Maxage) # Resident Leslie
matrix
Invader.matrix <- LESLIE(X.invader, Maxage) # Invader leslie
matrix
F.resident <- Resident.matrix[1,] # Resident DI fertility
F.invader <- Invader.matrix[1,] # Invader DI fertility
Maxgen <- 30 # Nos of gens to run
n.resident <- matrix(0,Maxage,1) # Resident population vector
n.resident[1] <- 1 # Initial resident popn size
for ( Igen in 2:Maxgen) # Iterate over generations
{
N <- sum(n.resident) # Total popn size
# Get DD fertility for resident population at time Igen
Resident.matrix[1,] <- DD.FUNCTION(ALPHA, BETA, F.resident, N)
n.resident <- Resident.matrix%*%n.resident # Resident popn
} # End of first Igen loop
### Introduce invader ####
Maxgen <- 100 # Number of generations to run
# Pre-allocate space for storage of invader population numbers
Pop.invader <- matrix(0,Maxgen,1)
# Pre-allocate space for invader vector
n.invader <- matrix(0,Maxage,1)
n.invader[1] <- 1 # Initial number of invaders
Pop.invader[1,1] <- n.invader[1] # Store initial numbers of
invaders
for (Igen in 2:Maxgen) # Iterate over generations
{
# Total number in population.
N <- sum(n.resident) þ sum(n.invader)
# DD fertility of resident
Resident.matrix[1,] <- DD.FUNCTION(ALPHA, BETA, F.resident, N)
# New resident vector
n.resident <- Resident.matrix%*%n.resident
# DD fertility of invader
Invader.matrix[1,] <- DD.FUNCTION(ALPHA, BETA, F.invader, N)
# New invader vector
n.invader <- Invader.matrix%*%n.invader
Pop.invader[Igen] <- sum(n.invader) # Store invader popn size
} # End of second Igen loop
# Now do linear regression of log(Pop.invader) on Generation
Generation <- seq(from¼1, to¼Maxgen) # Generate generation vector
Nstart <- 20 # Generations to ignore
# Linear regression
Invasion.model <-
lm(log(Pop.invader[Nstart: Maxgen])Gener ation[Nstart:Max gen])
INVASIBILITY ANALYSIS 191
# Elasticity value ¼ regression slope
Elasticity <- Invasion.model$coeff[2]
} # End of POP.DYNAMICS function
####################### MAIN PROGRAM #######################
N1 <- 30 # Nos of increments
X.Resident <- seq(from¼1, to¼3, length¼N1) # Resident body sizes
X.Invader <- X.Resident # Invader body sizes
d <- expand.grid(X.Resident, X.Invader) # Combinations
# Generate values at combinations
z <- apply(d,1,POP.DYNAMICS)
z.matrix <- matrix(z, N1, N1) # Convert to a matrix
# Plot contours
contour(X.Resident, X.Invader,z.matrix, xlab¼“Resident”,
ylab¼“Invader”)
# Place circle at predicted optimal body size
points( 1.68703, 1.68703, cex¼3) # cex triples size of circle
OUTPUT: (Figure 3.5)
Resident
Invader
1.0
1.0
1.5
2.0
2.5
3.0
1.5 2.0 2.5 3.0
–0.14
–0.12
–0.1
–0.08
–0.06
–0.06
–0.04
–0.04
–0.02
–0.08
0.06
0.04
0.02
0.02
0.04
0.06
0.08
0.01
0
0
0
0
–0.02
Figure 3.5 Pairwise invasibility plot for Scenario 2. The circle indicates the value obtained
from the elasticity analysis.
192 MODELING EVOLUTION
There is a single putative ESS, which, as shown by the circle, is confirmed by the
elasticity analysis described below.
3.4.5 Elasticity analysis
The program consists of the same components with the following changes
1. POP.DYNAMICS(X, Coeff): The body size of the resident and the multiplier
for the invader is passed to POP.DYNAMICS. As suggested by Benton and Grant
(1999), the value of the invader trait is set at 0.995 times that of the resident
(Coeff ¼ 0.995 in call to POP.DYNAMICS and hence X.invader ¼ 0.995*X.
resident). The change in population size of the invader population, Pop.
invader is stored and after the specified number of generations (Maxgen)
the elasticity is estimated as the slope of the linear regression of log(Pop.
invader)on Generation. To avoid poor estimates due to the initial stabiliza-
tion of the invader population the regression ignores the first 20 generations.
2. The optimum body size is that value at which elasticity is zero. This is estimated
first by calling the R function uniroot and then visually checked by plotting
elasticity versus body size and also the invasion exponent versus body size.
R CODE:
rm(list¼ls()) # Clear workspace
LESLIE <- function(x,Maxage) # Function to generate Leslie
matrix
{
CODING SAME AS IN INVASIBILITY ANALYSIS
} # End of Leslie function
############### Density-dependence function ###############
DD.FUNCTION <- function(ALPHA,BETA,Fi,n) {Fi*ALPHA*exp
(-BETA*n)}
###############Population dynamics function ###############
POP.DYNAMICS <- function(X, Coeff)
{
X.Resident <- X # Body size of resident population
X.invader <- X.Resident*Coeff # Body size of invader
REST OF CODE SAME AS IN INVASIBILITY ANALYSIS
} # End of POP.DYNAMICS function
####################### MAIN PROGRAM #######################
par(mfrow¼c(2,2)) # Divide graphics page into quarters
# Plot elasticity vs x
N.int <- 20 # Number of increments
X <- matrix(seq(from¼.5, to¼3, length¼N.int)) # Sequence of
body sizes
# Calculate elasticities for sequence
Elasticity <- apply(X,1,POP.DYNAMICS,0.995)
INVASIBILITY ANALYSIS 193
plot(X, Elasticity, type¼’l’) # Plot elasticity as a function of X
lines(c(.5,3), c(0,0)) # Add horizontal line at zero
# Calculate the optimum by calling uniroot
Optimum <- uniroot(POP.DYNAMICS, interval¼c(0.5,3),0.995)
Best.X <- Optimum$root # Save optimum X
print(c(“Optimum body size ¼”,signif(Best.X,6))) # Print out value
# Now plot Invasion exponent when resident is optimal
# Note that because of order in call cannot use apply here
# Convert the X sequence to coefficients for call to POP.DYNAMICS
Coeff <- X/Best.X
Invasion.exponent <- matrix(0,N.int,1) # Pre-allocate space
# Loop through values of X comparing to Best.X
for (i in 1:N.int){Invasion.exponent[i] <- POP.DYNAMICS(Best.X,
Coeff[i]) }
plot(X, Invasion.exponent, type¼’l’) # Plot invasion exponent vs x
points(Best.X, 0, cex¼2) # Plot point at predicted
optimum
OUTPUT: (Figure 3.6)
[1] “Optimum body size ¼” “1.68703”
The results from uniroot indicate that the optimum under density-dependent
regulation is smaller than in the density-independent case (Scenario 1). The plot
of invasion exponent on body size confirms that the putative ESS is indeed
an ESS.
3.5 Scenario 3: Functional dependence in the Ricker model
Ebenman et al. (1996) studied a stage-structured model in which selection favors
stability, whereas oscillatory behavior is favored in the age-structure model
XX
0.5
Elasticity
Invasion. exponent
–0.002
0.000
0.002
1.0 1.5 2.0 2.5 3.0
0.5
–0.30 –0.20 –0.10 0.00
1.0 1.5 2.0 2.5 3.0
Figure 3.6 Graphical analysis of elasticity and invasion exponent as a function of body
size as determined by an analysis of Scenario 2. The circle demarks the position of
the predicted optimum.
194 MODELING EVOLUTION
studied by Greenman et al. (2005). In this model we consider how evolution will
shape population dynamics by an analysis of the optimal parameter values in the
Ricker model. No age or stage structure is assumed.
3.5.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.
3. The parameters of the recruitment function are related such that the density-
independent component is negatively related to the density-dependent
component.
3.5.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
ð3:28Þ
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 ð3:29Þ
Examples of the population dynamics for increasing values of a are shown in
Figure 3.7. For low values of a the population reaches a stable equilibrium but as a
is increased the dynamics first become cyclical and then chaotic.
3.5.3 Pairwise invasibility analysis
The coding follows the general pattern of that in Scenario 2, again plotting the
value from the subsequent elasticity analysis on the contour plot. The DD.FUNC-
TION passes back the new population size using equation (3.28).
R CODE:
rm(list¼ls()) # Clear memory
DD.FUNCTION <- function(ALPHA,N1,N2) # Density-dependence
function
{
BETA <- ALPHA*0.001 # Set value of beta
N <- N1*ALPHA*exp(-BETA*N2) # New population size
return(N)
} # End of DD.FUNCTION
INVASIBILITY ANALYSIS 195
N(t)
N(t+1)
0 0 20 40 60 80 100
0 100 200 300
N(t+1)
0 100 200 300
N(t+1)
0 100 200 300
N
0 100 200 300
N
0 100 200 300
N
0 100 200 300
200 400 600 800 1,000
N(t)
0 200 400 600 800 1,000
Generation
0 20406080100
Generation
N(t)
0 200 400 600 800 1,000
0 20406080100
Generation
Figure 3.7 Population dynamics in the Ricker model in which b = 0.001a. From top to
bottom the values of a are 2, 10, and 20.
R CODING:
rm(list=ls()) # Clear workspace
par(mfrow=c(3,2)) # Divide graphics page into 3x2 panels
DD.FUNCTION <- function(n,ALPHA, BETA) {ALPHA*exp(-BETA*n)}
# main program
A <- c(2,10,20) # Values of alpha
for( j in 1:3) # Iterate over values of alpha
{
N.t <- seq(from=0, to=1000) # Population sizes
ALPHA <-A[j] # alpha
BETA <- ALPHA*0.001 # Beta
# Plot N(t+1) vs N(t)
N <- length(N.t) # Nos of values of N(t)
N.tplus1 <- matrix(0,N) # Pre-allocate space for N(t+1)
for ( i in 1:N) # Iterate over values of N
{
N.tplus1[i] <- N.t[i]*DD.FUNCTION(N.t[i],ALPHA,BETA)
} # End of N(t+1) on N(t) calculation
plot(N.t, N.tplus1, type=’l’, xlab=’N(t)’,ylab=’N(t+1)’)
lines(N.t, N.t) # Plot the line of equality
# Plot N(t) vs t
Maxgen <- 100 # Number of generations
N <- matrix(0, Maxgen) # Pre-allocate space for N(t)
N[1] <- 1 # Initial vale of N
for ( Igen in 2:Maxgen) # Iterate over generations
{
N[Igen] <- N[Igen-1]*DD.FUNCTION(N[Igen-1],ALPHA,BETA)
} # End of Igen loop
Generation <- seq(1,Maxgen) # Vector of generation numbers
plot(Generation, N, type=’l’) # Plot population trajectory
} # End of j loop
#### Function specifying population dynamics ####
POP.DYNAMICS <- function(ALPHA)
{
ALPHA.resident <- ALPHA[1] # Alpha for resident
ALPHA.invader <- ALPHA[2] # Alpha for invader
Maxgen1 <- 50 # Generations when only invader
present
Maxgen2 <- 300 # Generations after invader
introduced
Tot.Gen <- Maxgen1þMaxgen2 # Total number of generations
N.resident <- N.invader <- matrix(0,Tot.Gen) # Allocate space
N.resident[1] <- 1 # Initial number of resident
N.invader[Maxgen1]<- 1 # Initial number of invader
for (Igen in 2:Maxgen1) # Iterate over only resident
{
N.resident[Igen] <- DD.FUNCTION(ALPHA.resident, N.resident
[Igen-1], N.resident[Igen-1])
} # End of resident only period
# Now add invader
J <- Maxgen1þ1 # Staring generation of this period
for (Igen in J:Tot.Gen) # Iterate after introduction of invader
{
N.total <- N.resident[Igen-1]þ N.invader[Igen-1] # Total popn
size
# Resident population size
N.resident[Igen] <- DD.FUNCTION(ALPHA.resident, N.resident
[Igen-1],N.total)
# Invader population size
N.invader[Igen] <- DD.FUNCTION(ALPHA.invader, N.invader[Igen-
1],N.total)
} # End of invasion period
Generation <- seq(from¼1, to¼Tot.Gen) # Generation sequence
Nstart <-10þ Maxgen1 # Starting point for regression
# Regression model
Invasion.model <- lm(log(N.invader[Nstart:Tot.Gen])~ Genera-
tion[Nstart:Tot.Gen])
Elasticity <- Invasion.model$coeff[2] # Elasticity
return(Elasticity)
} # End of POP.DYNAMICS function
############# MAIN PROGRAM #############
N1 <- 30 # Nos of increments
A.Resident <- seq(from¼2, to¼ 4, length¼
N1) # Resident alpha
A.Invader <-
A.Resident # Invader alpha
INVASIBILITY ANALYSIS 197