Roff D.A. Modeling Evolution: An Introduction to Numerical Methods
Подождите немного. Документ загружается.
5. Up to three kills can be made per day: designate the number of kills as k, where
k ¼ 0, 1, 2, and 3.
6. The probability of each kill is independent and hence the number of kills per
day is a binomial variable:
PZ¼
kY
i
¼ Pkill
i;k
¼
3
k
p
k
i
ð1 p
i
Þ
3k
ð6:16Þ
where Z is the amount per individual per day.
7. An animal whose gut contents fall to or below Xcritical (i.e., x Xcritical) is dead.
8. Gut capacity has a maximum value of Xmax
.
9. Fitness is the survival probability to the end of some specified time interval, T.
6.5.3 Outcome chart and expected lifetime fitness function
To equate this scenario with the patch-foraging model, we note that pack size is
equivalent to patch identity and thus for each pack size, i, and kill number, k:
Benefit is a function of pack size and number of kills. We require two matrices, one
for the benefits and one for the probability, in which the rows correspond to pack
size and the columns to kills. For simplicity, we shall set pack size from 1 to 4. The
probabilities of a kill in relation to these pack sizes are 0.15, 0.31, 0.33, and 0.33,
respectively. Prey size is set at 11.25 kg.
R CODE:
rm(list¼ls()) # Remove all objects from memory
Benefit <- matrix(0,4,4) # Rows ¼ pack size, Columns ¼ number
of killsþ1
Pbenefit <- matrix(0,4,4) # Rows ¼ pack size, Columns ¼ number
of killsþ1
Pi <- c(0.15, 0.31, 0.33, 0.33) # Probability of single
kill for pack size
Y <- 11.25 # Size of single prey
k <- c(0,1,2,3) # Number of kills
for ( PackSize in 1:4) # Iterate over pack sizes
{
# Calculate binomial probabilities using function dbinom
Pbenefit[PackSize,] <- dbinom(x¼k, size¼3, prob¼Pi[Pack-
Size])
Kill made x > Xcritical Survives X
Yes — Yes min(x – Cost + Benefit
i,k
, Xmax)
............................................................................................................
No Yes Yes x – Cost
............................................................................................................
No Yes No Set at Xcritical
............................................................................................................
No No — Set at Xcritical
368 MODELING EVOLUTION
# Calculate benefits
Benefit[PackSize, 2:4] <- k[2:4]*Y/PackSize
}
Benefit # Print out Benefit matrix
Pbenefit # Print out Probability matrix
OUTPUT: (Table 6.8)
The MATLAB coding is given in the coding for the entire program (see section 6.9.5).
The value of Xcritical is assumed to be zero, that is, animals with no gut contents
have starved to death. This presents a minor difficulty in coding as now the index
value does not correspond to the state value. In this case, given that Xmax ¼ 30, we
could simply raise Xcritical to 1. However, to illustrate the general approach we
shall here retain Xcritical ¼ 0 and use the method of index adjustment previously
given. The parameter Xmin is replaced by Xinc, which is used to translate from x
to the index value. Note that the maximum value of the index, Max.Index,is
passed to function OVER.STATES.
A second complication is that changes in x do not follow unit steps: I shall follow
the suggestion of Mangel and Clark (1988) and use linear interpolation as dis-
cussed above and graphically illustrated in Figure 6.2. The coding for this is given
in the function FITNESS. The parameter Pmortality has been deleted as it is not
used.
Survival (¼ fitness) for each pack size is given by
Wðx; i; t 1Þ¼
X
3
k¼0
Pbenefit
i;k
Fðx; tÞð6:17Þ
and the optimal pack size is that pack size which maximizes survival, the fitness
being
Fðx; t 1Þ¼maxfWðx; 0; tÞ; Wðx; 1; tÞ; ...; Wðx; 4; tÞg ð6:18Þ
The final complication is that there may be several pack sizes that give the same
fitness. I shall not here consider whether there is a biologically reasonable way to
resolve this question but deal with the problem of locating those transitions
Table 6.8 Matrices showing benefits and probability of benefits as a function of pack size
Pack Size
Number of kills + 1 Probablity
[,1] [,2] [,3] [,4] [,1] [,2] [,3] [,4]
[1,] 0.00 11.25 22.5 33.75 0.6141 0.3251 0.0574 0.0034
......................................................................................................................................
[2,] 0.00 5.625 11.25 16.875 0.3285 0.4428 0.1989 0.0298
......................................................................................................................................
[3,] 0.00 3.75 7.5 11.25 0.3008 0.4444 0.2189 0.0359
......................................................................................................................................
[4,] 0.00 2.8125 5.625 8.4375 0.3008 0.4444 0.2189 0.0359
Note: R output is slightly modified for clarity.
DYNAMIC PROGRAMMING 369
where this occurs. At the start of the simulation we create a matrix called
CHOICES that holds a flag indicating whether there are several equivalent choices:
arbitrarily I designate 0 to indicate only one choice and 1 to indicate more than
one equivalent choice. At the commencement of the simulation all cells are set to
0. The actual test for equivalent choices is done in the function OVER.PATCHES.
The matrix passed back from this function has two columns and so we create a
1 2 vector called Choice that consists of two zeros:
Choice <- c(0,0)
one of these simply being a dummy to permit concatenation. Next we test if the
fitnesses in the first two rows of the sorted values are the same and if they are the
same we set the values in Choice to ones:
if(RHS[Best[1]]¼¼RHS[Best[2]]) Choice <- c(1,1)
This vector is then added to the bottom of the matrix Temp that is passed back to
OVER.STATES:
Temp <- rbind(Temp,Choice)
The data are then extracted, passed back to the main program, and stored in
CHOICES.
6.5.4 Calculating the decision matrix
R CODE:
rm(list¼ls()) # Remove all objects from memory
# Function to calculate fitness when organism is in state X
FITNESS <- function(X, Xcritical, Xmax, Xinc, Cost, Benefit,
Pbenefit, F.vectors)
{
# Note that the state value X is passed
# Note also that in this function Benefit and Pbenefit are vectors
# Iterate over the four kill values (0,1,2,3)
Max.Index <-1þ (Xmax-Xcritical)/Xinc # Get maximum index value
W <-0 # Set Fitness to zero
Xstore <- X # Set X to Xstore to preserve value through loop
for (I.Kill in 1:4) # Begin loop
{
X <- Xstore - Cost þ Benefit[I.Kill] # Calculate new state
value
# If X greater than Xmax then X must be set to Xmax
X <- min(X, Xmax)
# If X less than or equal to Xcritical then set to Xcritical
X <- max(X, Xcritical)
# Convert to Index value
370 MODELING EVOLUTION
Index <-1þ(X-Xcritical)/Xinc
# Index value probably not an integer
# So consider two integer values on either size of X
Index.lower <- floor(Index) # Choose lower integer
Index.upper <- Index.lower þ 1 # Upper integer
# Must stop index exceeding Max.Index. Not that Qx¼0 in this case
Index.upper <- min(Index.upper, Max.Index)
Qx <- X – floor(X) # qx for linear interpolation
W <-Wþ Pbenefit[I.Kill]*(Qx*F.vectors[Index.upper,2]þ(1-Qx)
*F.vectors[Index.lower,2])
} # End of I.Kill loop
return(W) # Return Fitness
} # End of function
# ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ....
# Function to iterate over patches i.e. over PACKS
OVER.PATCHES <- function(X, F.vectors, Xcritical,Xmax, Xinc,
Npatch, Cost, Benefit, Pbenefit)
{
RHS <- matrix(0,Npatch,1) # Set matrix for Right Hand
Side of equn
for (i in 1: Npatch) # Cycle over patches ¼ pack
sizes
{
# Call Fitness function. Pass Benefit and Pbenefit as vectors
RHS[i] <- FITNESS(X, Xcritical, Xmax, Xinc, Cost, Benefit[i,],
Pbenefit[i,], F.vectors)
} # End of i loop
# Now find optimal patch Best row is in Best[1]
Best <- order(RHS, na.last¼TRUE, decreasing¼TRUE)
Index <-1þ(X-Xcritical)/Xinc # Get Index value
F.vectors[Index,1]<- RHS[Best[1]] # Get best W ¼ F(x,t,T)
# Get best patch (¼pack). Remember to convert from index value
Best.Patch <- Best[1]
# Concatenate F(x,t) and the optimal patch (¼pack) number
Temp <- c(F.vectors[Index,1], Best.Patch)
# Add Temp to bottom of F.vectors and rename to Temp
Temp <- rbind(F.vectors, Temp)
# Create 1x2 vector to hold decision on more than one choice
# We only need one cell but it is convenient to use 2 for concatena-
tion onto Temp, as indicated below
# Set Choice to zero
Choice <- c(0,0)
if(RHS[Best[1]]¼¼
RHS[Best[2]]) Choice <-
c(1,1) # Equal
fitnesses
DYNAMIC PROGRAMMING 371
Temp <- rbind(Temp,Choice) # Bind to bottom of matrix
return (Temp)
} # End of function
# ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ....
# Function to iterate over states of X
OVER.STATES <- function(F.vectors, Xcritical, Xmax, Xinc, Npatch,
Cost, Benefit, Pbenefit, Max.Index)
{
Store <- matrix(0,Max.Index,3) # Create matrix for output
for ( Index in 2 : Max.Index) # Iterate over states of X
{
# For given X call Over.Patches to determine F(x,t) and best patch
X <- (Index-1)*Xinc þ Xcritical
Temp <- OVER.PATCHES(X, F.vectors, Xcritical, Xmax, Xinc,
Npatch, Cost, Benefit, Pbenefit)
# Extract components. Penultimate row is F(x,t,T) and best patch
n <- nrow(Temp)-2
F.vectors <- Temp[1:n,]
Store[Index,1:2] <- Temp[nþ1,] # Save F(x,t,T) and best patch
Store[Index,3] <- Temp[nþ2,1] # Save Flag for several choices
} # End of X loop
# Add Store values to end of F.vectors for pass back to main program
Temp <- cbind(F.vectors, Store) # Combined by columns
return(Temp) # Return F.vectors and Store
} # End of function
# ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ....
# MAIN PROGRAM
# Initialize parameters
Xmax <- 30 # Maximum value of X ¼ gut capacity
Xcritical <- 0 # Value of X at which death occurs
Xinc <- 1 # Increment in state variable
Max.Index <-1þ (Xmax-Xcritical)/Xinc # Maximum index value
Cost <- 6 # Cost ¼ Daily food requirement
Npatch <- 4 # Number of patches¼ packs
# Calculate benefit as a function of pack size (rows)
# and number of kills (columns)
Benefit <- matrix(0,4,4) # Rows ¼ pack size, Columns
¼ number of killsþ1
Pbenefit <- matrix(0,4,4) # Rows ¼ pack size, Columns ¼
number of killsþ
1
#
Probability of single kill for pack size
Pi <- c(0.15, 0.31, 0.33, 0.33)
Y <- 11.25 # Size of single prey
372 MODELING EVOLUTION
k <- c(0,1,2,3) # Number of kills
for ( PackSize in 1:4) # Iterate over pack sizes
{
# Calculate binomial probabilities using function dbinom
Pbenefit[PackSize,] <- dbinom(x¼k, size¼3, prob¼Pi
[PackSize])
# Calculate benefits ¼ amount per individual
Benefit[PackSize, 2:4] <- Y*k[2:4]/PackSize
}
Horizon <- 31 # Number of time steps
# Set up matrix for fitnesses
# Column 1 is F(x, t). Column 2 is F(x, tþ1)
F.vectors <- matrix(0, Max.Index,2)
F.vectors[2:Max.Index,2] <- 1 # Cell 1,2 ¼ 0 ¼ Dead
# Create matrices for output
FxtT <- matrix(0,Horizon,Max.Index) # F(x,t,T)
Best.Patch <- matrix(0,Horizon,Max.Index) # Best patch
number
# Matrix for flag indicating multiple equivalent choices
#0¼ only one choice, 1 ¼ more than one choice
CHOICES <- matrix(0,Horizon,Max.Index)
# Start iterations
Time <- Horizon # Initialize Time
while ( Time > 1)
{
Time <- Time - 1 # Decrement Time by 1 unit
# Call OVER.STATES to get best values for this time step
Temp <- OVER.STATES(F.vectors, Xcritical, Xmax, Xinc,
Npatch, Cost, Benefit, Pbenefit, Max.Index)
# Extract F.vectors
TempF <- Temp[,1:2]
# Update F1
for ( J in 2: Max.Index) { F.vectors[J,2] <- TempF[J,1]}
# Store results
Best.Patch[Time,] <- Temp[,4]
FxtT[Time,] <- Temp[,3]
CHOICES[Time,] <- Temp[,5]
} # End of Time loop
# Output information. For display add states to last row of matrices
# Note that state variable conversion from index value
Index <- seq(from¼1, to¼Max.Index)
Best.Patch[Horizon,] <-
(Index-1)*XincþXcritical
FxtT[Horizon,] <- (Index-1)*XincþXcritical
Best.Patch[,1:Max.Index] # Print Decision matrix
DYNAMIC PROGRAMMING 373
signif(FxtT[,1:Max.Index],3) # Print Fxt of Decision
matrix: 3 sig places
CHOICES[,1:Max.Index] # Print matrix indicating choice flag
# Plot data
y <- Best.Patch[Horizon,2:Max.Index]
x <- seq(from¼1, to ¼Horizon-1)
par(mfrow¼c(2,2))
persp(x, y, Best.Patch[1:30,2:Max.Index], xlab¼‘Time’,
ylab¼‘x ¼ Gut contents’, zlab¼‘Optimal Pack size’, theta¼ 20,
ph¼25, lwd¼ 1) # 3D plot
image(x, y, Best.Patch[1:30,2:Max.Index], col¼terrain.colors
(50), xlab¼‘Time’, ylab¼‘x ¼ Gut contents’, las¼1) # Colored grid
image(x, y, CHOICES[1:30,2:Max.Index], col¼terrain.colors
(50), xlab¼‘Time’, ylab¼‘x ¼ Gut contents’, las¼ 1) # Colored grid
OUTPUT: (Figure 6.3)
One Pack size optimal
Pack size = 1
Pack size = 2
5
5
10
15
x = Gut contents
x = Gut contents
20
25
30
5
10
15
x = Gut contents
20
25
30
Time
Optimal Pack size
10 15
Time
20 25 30
51015
Time
20 25 30
Figure 6.3 Results for Scenario 3: Top shows two graphical representation of the
decision matrix (unlabeled color = optimal pack size of 3) and bottom a visualization of
the choice matrix (unlabeled color = more than one optimal choice).
374 MODELING EVOLUTION
The decision matrix, the matrix of fitnesses, and the matrix indicating the pres-
ence of multiple equivalent choices are printed out but not shown here. Figure 6.3
shows two visualizations of the decision matrix and a visualization of the matrix
CHOICE. Over most of state space a pack size of 2 is optimal. The number of cases
in which there are multiple equivalent choices increase with the state value and
the approach of the end of the time span.
MATLAB CODE: see Section 6.9.5
6.6 Scenario 4: Host choice in parasitoids: fitness
decreases with time
A frequent use of dynamic programming is to examine oviposition behavior in
organisms such as parasitoids that lay clutch sizes that depend upon host or patch
quality. The important change in this scenario compared to the previous ones is that
thevalueofthestatevariableincreasesaswemovetowardt ¼ 1, rather than
decreasing.
6.6.1 General assumptions
1. The animal commences the time period with some fixed quantity of eggs, as
occurs, for example, in some Lepidoptera. In general, animals can be classified
into capital breeders that use only, or primarily, resources gathered prior to
maturity and income breeders that garner resources for reproduction after
maturity. The present model applies to capital breeders, though it can easily
be adapted for income breeders.
2. Patches or hosts vary in quality.
3. The survival and growth of larvae depend on the number in the clutch and host
quality.
4. Variation in host quality can be detected by the ovipositing females.
5. Survival of the female may or may not change over time. For computational
simplicity, we assume that the sequence of events is that egg-laying precedes
the determination of survival over the time period.
6. Only one host at most is encountered per time interval.
7. Hosts already with eggs are not encountered.
8. Fitness is a function of the number of offspring.
6.6.2 Mathematical assumptions
1. There are four types of host.
2. The single host fitness can be modeled by a cubic function:
DYNAMIC PROGRAMMING 375
Benefit
i;n
¼ a
i
þ b
i
n þ c
i
n
2
þ d
i
n
3
ð6:19Þ
where n is the number of eggs laid on a host, the subscript i refers to host type, and
the coefficients vary according to host type. In this case we have not explicitly
defined fitness in relation to the components of offspring survival and future
reproduction, but have absorbed these into a single function empirically derived
by Charnov and Skinner (1984) for the parasitoid wasp Nasonia vitipennis.An
important feature of this function is that it has an intermediate optimum, but
the fitness curve for host type 1 is clearly incorrect and derives from the fact that
the model is extended beyond the observed range (Figure 6.4). It will never be
optimal to increase clutch size beyond the local maximum (see below for
Cluth size
Female Fitness
0
0
5
10
15
20
25
30
10 20 30 40
Figure 6.4 Fitness increments on each type of host parasitized by the wasp Nasonia
vitripennis modeled by a cubic function. Coefficients from Table 4.1 of Mangel and Clark
(1988). Coding to generate Benefits matrix and plot data prior to setting values greater
than n* (the single host maximum) to zero. Because zero occupies the first column we
apply an index transformation.
rm(list=ls()) # Remove all objects from memory
Xmax <- 40 # Maximum value of X = eggs
Xcritical <- 0 # Lowest value of X = 0 eggs
Xinc <- 1 # Increment in state variable
Max.Index <- 1 + (Xmax-Xcritical)/Xinc # Max Index value
# Create host coefficient matrix from which to get Benefits
Host.coeff <- matrix(0,4,4)
Host.coeff[1,] <- c(-0.2302, 2.7021, -0.2044, 0.0039)
Host.coeff[2,] <- c(-0.1444, 2.2997, -0.1170, 0.0013)
376 MODELING EVOLUTION
derivation) and hence all values greater than this can be set to zero (coding given
in figure caption). A plausible model for this type of function (i.e., single maxi-
mum) is that offspring survival and body size, which controls future fecundity,
decreases with clutch size but, because fitness is equal to clutch size times, the
expected fecundity of each offspring, fitness initially increases with clutch size.
The optimum clutch size for a single clutch can be obtained from the calculus
dBenefit
i;n
dn
¼ b
i
þ 2c
i
n þ 3d
i
n
2
ð6:20Þ
The optimum clutch size is then found by setting dBenefit
i;n
=dn ¼ 0 and solving the
resultant quadratic (see Scenario 1 of Chapter 2), say n*. It will never be optimal for
a female to lay more eggs than n*, but it could be optimal to lay fewer eggs if the
host is of poor quality and the female is likely to find higher quality hosts in the
future. As noted above, to avoid the unreal behavior of at least one of the single
host fit we set values greater than n* equal to zero.
3. The probability of encountering a host is constant but different for each host,
designated as Pbenefit
i
, where i is the ith host type. Thus the probability of not
encountering a host, P
0
,is
Host.coeff[3,] <- c(-0.1048, 2.2097, -0.0878, 0.0004222)
Host.coeff[4,] <- c(-0.0524, 2.0394, -0.0339, -0.0003111)
# Calculate benefit as a function of
# clutch size (rows) and Host type (columns)
Clutch <- seq(from = 0, to = Xmax)
Benefit <- matrix(0, Max.Index, 4) # Zero to Xmax
for (I.Host in 1:4) # Iterate over host types
{
Benefit[,I.Host] <- Host.coeff[I.Host,1] + Host.coeff[I.
Host,2]*Clutch + Host.coeff[I.Host,3]*Clutch^2 + Host.coeff[I.
Host,4]*Clutch^3
}
# Plot data
plot(Clutch, Benefit[,1],type=‘l’, xlab=‘Clutch size’,ylab=“-
Female Fitness”, las=1, lwd=4)
# lwd = line width, lty = line type,1=solid, 2=dashed, 3=dotted,
4=dotdash,
for (i in 2:4){lines(Clutch, Benefit[,i],type=‘l’, lwd=4,
lty=i)}
SHM <- c(9,12,14,23) # Set single host maximum. See text for
derivation
# Make all values > than SHM=0. Note that we use 1 because of zero
class
for (i in 1:4){Benefit[(SHM[i]+1):Max.Index,i] <- 0}
DYNAMIC PROGRAMMING 377