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

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

A <- c(100, 1) # a coefficients

Xegg <- c( 2, 1) # “optimum” egg size

Xclutch <- c(5, 10) # “Optimum” clutch size

P1 <- 0.5 # Probability of host 1

Psurvival <- 0.90 # Survival probability per time increment

Npatch <- 2 # Number of patches¼ hosts

Horizon <- 10 # Number of time steps

# Set up matrix for fitnesses

# Column 1 is temporary F(x, tþ1). Column 2 is F(x, tþ1) Both are

zero

F.vectors <- matrix(0, Max.Index,2)

# Create matrices for output

FxtT <- matrix(0,Horizon,Max.Index)# F(x,t)

Best.Clutch1 <- matrix(0,Horizon,Max.Index)# Best clutch size

for host 1

Best.Clutch2 <- matrix(0,Horizon,Max.Index)# Best clutch size

for host 2

Best.Egg1 <- matrix(0,Horizon,Max.Index)# Best egg size for host 1

Best.Egg2 <- matrix(0,Horizon,Max.Index)# Best egg size for host 2

Choice.H1 <- matrix(1,Horizon,Max.Index)# 1 or 2 choices for

host type 1

Choice.H2 <- matrix(1,Horizon,Max.Index)# 1 or 2 choices for

host type 2

# 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,

Psurvival, Max.Index, Wmax, A, Xegg, Xclutch, P1)

# Extract F.vectors

TempF <- Temp[,1:2]

# Update F1

for ( J in 2: Max.Index) { F.vectors[J,2] <- TempF[J,1]}

# Store results

FxtT[Time,] <- Temp[,3]

Best.Clutch1[Time,] <- Temp[,4]

Best.Egg1[Time,] <- Temp[,5]

Best.Clutch2[Time,] <- Temp[,6]

Best.Egg2[Time,] <- Temp[,7]

Choice.H1[Time,] <

- Temp[,8]

Choice.H2[Time,] <-

Temp[,9]

} # End of Time loop

# Output information. For display add states to last row of matrices

398 MODELING EVOLUTION

Index <- seq(from¼1, to¼Max.Index)

FxtT[Horizon,] <- (Index-1)*XincþXcritical

Best.Clutch1[Horizon,] <- (Index-1)*XincþXcritical

Best.Clutch2[Horizon,] <- (Index-1)*XincþXcritical

Best.Egg1[Horizon,] <- (Index-1)*XincþXcritical

Best.Egg2[Horizon,] <- (Index-1)*XincþXcritical

Choice.H1[Horizon,] <- (Index-1)*XincþXcritical

Choice.H2[Horizon,] <- (Index-1)*XincþXcritical

# Print Best clutch, best egg choice flag and decision matrix

Best.Clutch1[,1:Max.Index] # Best clutch on host 1

Best.Clutch2[,1:Max.Index] # Best clutch on host 2

signif(Best.Egg1[,1:Max.Index],3) # Best egg on host 1

signif(Best.Egg2[,1:Max.Index],3) # Best egg on host 2

Choice.H1[,1:Max.Index] # 1¼only 1 choice 2 ¼1þ

Choice.H2[,1:Max.Index] # 1¼only 1 choice 2 ¼1þ

signif(FxtT[,1:Max.Index],3) # Print Fxt of Decision

matrix: 3 sig places

OUTPUT:

The optimum egg–clutch size combinations are shown in Table 6.9. Multiple

optima occur only at the very lowest value of reproductive biomass. Variation in

clutch size is evident across reproductive biomass but not over the time period.

Egg size varies primarily with reproductive biomass but not time. When repro-

ductive biomass is relatively low the optimum egg size tends to be larger than the

single host optimum (cf. Figure 6.9).

MATLAB CODE: See Section 6.9.8.

6.8 Some exemplary papers

Harris, W. E. and J. R. Lucas. 2002. A state-based model of sperm allocation in

a group-breeding salamander. Behavioral Ecology 13:705–712.

Question: How should males of the small-mouthed salamander allocate spermato-

phores in the face of competitors and a variable number of females?

State variable: Number of spermatophores a male has available at any time during

the breeding season. Spermatophore number is set at the start of the breeding

season.

Fitness: Maximization of reproductive success ¼ Number of clutches fertilized

Terminal ﬁtness: F(x, T) ¼ 0

Forward simulation: Yes

Experimental test: Yes

DYNAMIC PROGRAMMING 399

Table 6.9

Optimal clutch and egg sizes as a function of reproductive biomass and time (bold

italic) on the two types of host modeled in Scenario 5

X=Time 123456789101

23 4 5 6 7 8 9 10

Optimum clutch size on host 1

Optimum clutch size on host 2

1112233444511

23456678

............................................................................................................................

..................................................................................................

12233444511234

56678

............................................................................................................................

..................................................................................................

122334445

1123456678

............................................................................................................................

..................................................................................................

12233444511

23456678

............................................................................................................................

..................................................................................................

5112233444511

23456678

............................................................................................................................

..................................................................................................

12233444511234

56678

............................................................................................................................

..................................................................................................

12223344511234

56678

............................................................................................................................

..................................................................................................

11223344511234

56678

............................................................................................................................

..................................................................................................

11223344511234

56678

............................................................................................................................

..................................................................................................

Optimum egg size on host 1

Optimum egg size on host 2

............................................................................................................................

..................................................................................................

2 1.5 2 1.67 2 1.75 2 2.04 2 3 2

1.5 1.33 1.25 1.2 1.17 1.33 1.26 1.25

............................................................................................................................

..................................................................................................

232 1.5 2 1.67 2 1.75 2 2.04 2 3 2

1.5 1.33 1.25 1.2 1.17 1.33 1.26 1.25

............................................................................................................................

..................................................................................................

2 1.5 2 1.67 2 1.75 2 2.04 2

3 2 1.5 1.33 1.25 1.2 1.17 1.33 1.26 1.25

............................................................................................................................

..................................................................................................

2 1.5 2 1.67 2 1.75 2 2.04 2 3 2

1.5 1.33 1.25 1.2 1.17 1.33 1.25 1.25

............................................................................................................................

..................................................................................................

532 1.5 2 1.67 2 1.75 2 2.04 2 3 2

1.5 1.33 1.25 1.2 1.17 1.33 1.24 1.25

............................................................................................................................

..................................................................................................

2 1.5 2 1.67 2 1.75 2 2.04 2 3 2 1.5 1.33

1.24 1.2 1.17 1.33 1.22 1.25

............................................................................................................................

..................................................................................................

2 1.5 2 2.05 2 2.04 2 2.03 2 3 2 1.45

1.33 1.2 1.2 1.17 1.24 1.16 1.25

............................................................................................................................

..................................................................................................

2 2.02 2 2.03 2 2.03 2 2.02 2 3 2 1.07

1.2 1.11 1.2 1.11 1.11 1.1 1.21

............................................................................................................................

..................................................................................................

22222222231

11111111

Note: Combinations in which simulation gave more than

one optimal choice are shown in bold font.

Peterson, J. H., B. D. Roitberg, and R. C. Ydenberg. 2007. When nesting in-

volves two sequential, mutually exclusive activities: what’s a mother to do?.

Evolutionary Ecology Research 9:1187–1197.

Question: When should hymenopteran parents stop provisioning the current nest

and decide whether to seal the entrance to the nest?

State variables:(a) Current pollen value collected and (b) current nest value

Fitness: Number of copies of alleles passed on (sons are haploid, females diploid).

Terminal ﬁtness: F(x

, x

, T) ¼ 0

Forward simulation: No

Experimental test :No

Kerkhoff, A. J. 2004. Expectation, explanation and masting. Evolutionary Ecolo-

gy Research 6:1003–1020.

Question: How should trees allocate energy to growth and reproduction?

State variables:(a) Mass and (b) stored reserves

Fitness: Total lifetime reproductive output

Terminal ﬁtness: F(x

, x

, T) ¼ 0

Forward simulation: Yes

Experimental test: No

Lessells, C. M. 2005. Why are males bad for females? Models for the evolution

of damaging male mating behavior. American Naturalist 165:S46–S63.

Question: Given that a male damages a female during mating, which reduces her

survival, what is the optimal female age-speciﬁc behavior (mate or do not mate;

oviposition rate)?

State variable: Level of damage

Fitness: Lifetime egg production

Terminal ﬁtness: F(x, T) ¼ 0

Forward simulation: Yes

Experimental test: No

Weber, T. P., B. J. Ens, and A. I. Houston. 1998. Optimal avian migration: A

dynamic model of fuel stores and site use. Evolutionary Ecology 12:377–401.

Question: What is the optimal decision (stay or depart) for a migrating bird at

location i at time t?

State variable: Amount of fuel as a percentage of lean body mass

Fitness: Function of amount of fuel at ﬁnal arrival, R(x), arrival time K(t), and a

constant, B

Terminal ﬁtness: R(x)K(T) þ B

Forward simulation: Yes

Experimental test: No

DYNAMIC PROGRAMMING 401

Galvani, A. and R. Johnstone. 1998. Sperm allocation in an uncertain world.

Behavioral Ecology and Sociobiology 44:161–168.

Question: What is the optimal allocation of sperm when females of varying quality

are encountered at random?

State variable:(a) Quantity of sperm and (b) female quality

Fitness: Cumulative payoff function that for each allocation decision depends on

sperm quantity and female quality

Terminal ﬁtness: F(x

, x

, T) ¼ 0

Forward simulation: Yes

Experimental test: No

Wajnberg, E., P. Bernhard, F. Hamelin, and G. Boivin. 2006. Optimal patch

time allocation for time-limited foragers. Behavioral Ecology and Sociobiology

60:1–10.

Question: What is the relationship between patch residence time and age for a

parasitoid searching for host on which to oviposit?

State variable: Age

Fitness: Accumulated number of offspring

Terminal ﬁtness: F(x, T) ¼ 0

Forward simulation: Yes

Experimental test: Yes

6.9 MATLAB code

6.9.1 An algorithm for constructing the decision matrix

Function to calculate fitness when organism is in state X

function W ¼FITNESS(X, Xcritical, Xmax, Xmin, Cost, Benefit, Pbe-

nefit, Pmortality, F_vectors)

% State in patch if forager finds food

X_Food ¼ X - Cost þ Benefit;

% If X_Food greater than Xmax then X_Food must be set to Xmax

X_Food ¼ min(X_Food, Xmax);

% If X_Food less than or equal to Xcritical then set to Xcritical

X_Food ¼ max( X_Food, Xcritical);

% State in patch if forager does not find food

X_NoFood ¼ X - Cost;

% If X_NoFood is less than Xcritical set X_NoFood to Xcritical

X_NoFood ¼ max(X_NoFood, Xcritical);

Term1 ¼ Pbenefit*F_vectors(X_Food,2); % If food is found

Term2 ¼(1-Pbenefit)*F_vectors(X_NoFood,2); % If food is not found

W ¼ (1 - Pmortality)*(Term1 þ Term2); % Survival in patch

% End of function

402 MODELING EVOLUTION

% ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

% Function to iterate over patches

function Temp ¼ OVER_PATCHES(X, F_vectors, Xcritical,Xmax, Xmin,

Npatch, Cost, Benefit, Pbenefit, Pmortality)

RHS ¼ zeros(Npatch); % Pre-allocate matrix for Right Hand Side of

equn

for i ¼ 1: Npatch; % Cycle over patches

% Call Fitness function

RHS(i) ¼ FITNESS(X, Xcritical, Xmax, Xmin, Cost, Benefit(i),

Pbenefit(i), Pmortality(i), F_vectors);

end % End of i loop

% Now find optimal patch C¼Highest RHS, I¼Row¼patch number

[C,I] ¼ max(RHS);

F_vectors(X,1) ¼ C(1);

Best_Patch ¼ I(1);

% Concatenate F(x,t,T) and the optimal patch number

Temp ¼ [F_vectors(X,1), Best_Patch];

% Add Temp to bottom of F.vectors and rename to Temp

Temp ¼ vertcat(F_vectors, Temp);

% End of function

% ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ....

% Function to iterate over states of X

function Temp¼OVER_STATES(F_vectors, Xcritical, Xmax, Xmin,

Npatch, Cost, Benefit, Pbenefit, Pmortality)

Store ¼ zeros(Xmax,2); % Create matrix for output

for X ¼ Xmin : Xmax % Iterate over states of X

% For given X call Over_Patches to determine F(x,t) and best patch

Temp ¼ OVER_PATCHES(X, F_vectors, Xcritical, Xmax, Xmin,

Npatch, Cost, Benefit, Pbenefit, Pmortality);

% Extract components_ Last row is F(x,t) and best patch

n ¼ size(Temp,1)-1;

F_vectors ¼ Temp(1:n,1:2);

Store(X,1:2) ¼ Temp(nþ1,1:2); % Save F(x,t,T) and best patch

end % End of X loop

% Add Store values to end of F_vectors for pass back to main program

Temp ¼ horzcat(F_vectors, Store); % Combined by columns

% End of function

% ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ....

% MAIN PROGRAM

% Initialize parameters

Xmax ¼ 10 ; % Maximum value of X

Xcritical ¼ 3; % Value of X at which death occurs

Xmin ¼ Xcriticalþ1; % Smallest value of X allowed

Cost ¼ 1; % Cost per period

Pmortality ¼ [0,

0.004, 0.02]; % Probability of mortality

DYNAMIC PROGRAMMING 403

Pbenefit ¼ [1, 0.4, 0.6]; % Probability of finding food

Benefit ¼ [0, 3, 5]; % Benefit if food is discovered

Npatch ¼ 3; % Number of patches

Horizon ¼ 20; % Number of time steps

% Set up matrix for fitnesses

% Column 1 is F(x, tþ1)_ Column 2 is F(x, t)

F_vectors ¼ zeros(Xmax,2);

F_vectors(Xmin:Xmax,2) ¼ 1;

% Create matrices for output

FxtT ¼ zeros(Horizon,Xmax); % F(x,t,T)

Best_Patch ¼ zeros(Horizon,Xmax); % Best patch number

% 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, Xmin, Npatch,

Cost, Benefit, Pbenefit, Pmortality);

% Extract F_vectors

TempF ¼ Temp(:,1:2);

% Update F1

for J ¼ Xmin: Xmax

F_vectors(J,2) ¼ TempF(J,1);

end % End of J loop

% Store results

Best_Patch(Time,:) ¼ Temp(:,4);

FxtT(Time,:) ¼ Temp(:,3);

end % End of Time loop

% Output information_ For display add wts to last row of matrices

X ¼ 1: Xmax;

Best_Patch(Horizon,:) ¼ X;

FxtT(Horizon,:) ¼ X;

Best_Patch(:,Xmin:Xmax) % Print Decision matrix

vpa(FxtT(:,Xmin:Xmax),3) % Print Fxt of Decision matrix: 3 sig

places

6.9.2 Using the decision matrix: individual prediction

Because the random number generator is not the same, the output will not be

exactly the same as for the R output. As with the R coding the decision matrix is

ﬁrst generated by running the previous program and then the following (do not

issue a clear all command):

% Initialize parameters

rand(‘twister’,10); % Set random number seed

404 MODELING EVOLUTION

Xmax ¼ 10; % Maximum value of X

Xcritical ¼ 3; % Value of X at which death occurs

Xmin ¼ Xcriticalþ1; % Smallest value of X allowed

Cost ¼ 1; % Cost per period

Pmortality ¼ [0, 0.004, 0.0]; % Probability of mortality

Pbenefit ¼ [1, 0.4, 0.6]; % Probability of finding food

Benefit ¼ [0, 3, 5]; % Benefit if food is discovered

Npatch ¼ 3; % Number of patches

Horizon ¼ 15; % Number of time steps

Output ¼ zeros(Horizon,10);% Matrix to hold output

Time ¼ 1: Horizon; % Values for x axis in plot

for Replicate ¼ 1: 10; % Iterate over 10 replicates

X ¼ 4; % Animal starts in state 4

for i ¼ 1:Horizon; % Iterate over time

if(X > Xcritical); % Check that animal still alive

Patch ¼ Best_Patch(i,X); % Select patch

% Check if animal survives predation

% Generate random number

if rand(1) < Pmortality(Patch)

‘Dead from predator at Replicate and time ¼ ’

[Replicate, i]

end % end if

% Now find new weight

% Set multiplier to zero, which corresponds to no food found

Index ¼ 0;

if rand(1) < Pbenefit(Patch)

Index ¼ 1; % food is discovered

end % end if

X ¼ X - Cost þ Benefit(Patch)*Index;

% If X greater than Xmax then X must be set to Xmax

X ¼ min(X, Xmax);

% If X less than X then animal dies

if( X< Xmin)

‘Dead from starvation at Replicate and time ¼’

[Replicate, i]

end % end if

Output(i,Replicate) ¼ Patch; % Store data

end % End of if(X > Xcritical)

end % End of time loop

% subplot divides the current figure into rectangular panes that

are

% numbered row-wise. Each pane contains an axes object.

% Subsequent plots are output to the current pane.

subplot(5,2,Replicate); % Divide graph page into 5x2 panels

plot(Time,

Output(:,Replicate))

DYNAMIC PROGRAMMING 405

xlabel(‘Time’); ylabel(‘Patch selected’)

axis([0 15 0 4])

end % End of replicate loop

6.9.3 Using the decision matrix: expected state

% Initialize parameters

Xmax ¼ 10; % Maximum value of X

Xcritical ¼ 3; % Value of X at which death occurs

Xmin ¼ Xcriticalþ1; % Smallest value of X allowed

Cost ¼ 1; % Cost per period

Time ¼ 2; % Current state to be considered

Pmortality ¼ [0, 0.004, 0.02]; % Probability of mortality

Pbenefit ¼ [1, 0.4, 0.6]; % Probability of finding food

Benefit ¼ [0, 3, 5]; % Benefit if food is discovered

% Set transition density matrix to zero

Trans_density ¼ zeros(Xmax, Xmax);

% Step 1 Cycle over all values of z from Xmin to Xmax

for z ¼ Xmin : Xmax % Iterate over states

% Select the best patch from the Decision matrix at row Time

K ¼ Best_Patch(Time,z); % Decision matrix is called Best_Patch

% Calculate w(x,t|z)

% Found food and survives predator

x ¼ min(z - Cost þ Benefit(K), Xmax);

% Assign probability

Trans_density(z,x)¼ (1-Pmortality(K))*Pbenefit(K);

% Food not found

x ¼ z - Cost;

% State exceeds the critical value

if x > Xcritical

% Animal survives

Trans_density(z,x)¼ (1-Pmortality(K))*(1-Pbenefit(K));

% Animal does not survive

Trans_density(z,Xcritical)= Pmortality(K);

% State is less than critical

else

Trans_density(z,Xcritical)= Pmortality(K)þ(1-Pmortality(K))

*(1-Pbenefit(K)); % All one line

end % End of if

end % end of z loop

Trans_density % Write out matrix

406 MODELING EVOLUTION

6.9.4 Scenario 2: Calculating the decision matrix

% Function to calculate fitness when organism is in state X

function W ¼FITNESS(X, Xcritical, Xmax, Xmin, Cost, Benefit,

Pbenefit, Pmortality, F_vectors)

% State in patch if forager finds food

X_Food ¼ X þ Benefit; % Eliminate Cost

% If X_Food greater than Xmax then X_Food must be set to Xmax

X_Food ¼ min(X_Food, Xmax);

% If X_Food less than or equal to Xcritical then set to Xcritical

X_Food ¼ max( X_Food, Xcritical);

% State in patch if forager does not find food

X_NoFood ¼ X; % Eliminate Cost

% If X_NoFood is less than Xcritical set X_NoFood to Xcritical

X_NoFood ¼ max(X_NoFood, Xcritical);

Term1 ¼ Pbenefit*F_vectors(X_Food,2); % If food is found

Term2 ¼ (1-Pbenefit)*F_vectors(X_NoFood,2); % If food is not

found

W ¼ (1 - Pmortality)*(Term1 þ Term2); % Survival in patch

% End of function

% ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

% Function to iterate over patches

function Temp ¼ OVER_PATCHES(X, F_vectors, Xcritical,Xmax,

Xmin, Npatch, Cost, Benefit, Pbenefit, Pmortality)

RHS ¼ zeros(Npatch,1); % Pre-allocate matrix for Right Hand Side

of eqn

for i ¼ 1: Npatch; % Cycle over patches

% Call Fitness function

RHS(i) ¼ FITNESS(X, Xcritical, Xmax, Xmin, Cost, Benefit(i),

Pbenefit(i), Pmortality(i,X), F_vectors);

end % End of i loop

% Now find optimal patch C¼Highest RHS, I¼Row¼patch number

[C,I] ¼ max(RHS);

F_vectors(X,1) ¼ C;

Best_Patch ¼ I;

% Concatenate F(x,t) and the optimal patch number

Temp ¼ [F_vectors(X,1), Best_Patch];

% Add Temp to bottom of F.vectors and rename to Temp

Temp ¼ vertcat(F_vectors, Temp);

% End of function

% ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

% Function to iterate over states of X

function Temp¼OVER_STATES(F_vectors, Xcritical, Xmax, Xmin,

Npatch, Cost, Benefit, Pbenefit, Pmortality)

DYNAMIC PROGRAMMING 407