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

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

% Now calculate x3 and fitness

x3 ¼ (R-N*(x1þx2))/N ; % Determine x2

if (x3 < 0) % Check if x3 exists(>0)

W¼0;

else % Check if first clutch mass exceeds reserves

if(N*x1 > R) % Propagule too large

W ¼ 0;

else

W1 ¼ N*S1*Fmax*(1-exp(-a*x1)); % Calculate first fecundity

% Calculate size of propagules in 2nd clutch and see if reserves

exceeded

if(N*x2 > R)

W ¼ W1; %Propagules in 2nd clutch too large

else

W2 ¼ N*S2*Fmax*(1-exp(-a*x2)); % Calculate 2nd fecundity

% Calculate the size of Propagules in 3rd clutch

% Note that there must be reserves remaining at this stage

W3 ¼ N*S3*Fmax*(1-exp(-a*x3)); % Calculate 3rd fecundity

W ¼ W1þW2þW3; % Sum fitness components

end % End 3rd else

end % End 2nd else

end % End 1st else

Propagules ¼[x1,x2,x3]; % Store sizes

save PROPAGULE.txt Propagules –ASCII % Output sizes

W¼ -W; % Return -fitness

Finally the main program:

clear all; % Clear the workspace

% Locate optima

fminsearch(@FITNESS,1); % Call fminsearch with estimate at 1

load PROPAGULE.txt % Get optima from file

PROPAGULE % Print out results

OUTPUT:

PROPAGULE ¼ 2.7705 1.2290 0.0005

As with R, MATLAB gives more or less the same answer as before.

2.18.30 Scenario 12: Finding the maximum using a numerical approach

Using nlm (R) or fminsearch (MATLAB)

function W¼FITNESS(v) % Function to evaluate fitness

{ This function is the same as that used in plotting}

158 MODELING EVOLUTION

Main Program:

clear all; % Clear the workspace

% MAIN PROGRAM

% Call fminsearch passing fitness function with initial estimates

X ¼ fminsearch(@FITNESS, [0.1,0.1]);% Store estimates in X

% Calculate x3

R ¼ 400; N ¼ 100; % Parameter values

x3 ¼ (R-N*(X(1)þX(2)))/N; % x3

vpa([X(1),X(2), x3],6) % Output propagule sizes

OUTPUT:

ans ¼[ 2.77073, 1.22927, .228709e-7]

Previous results:

ans ¼ 2.7705 1.2291 0.0005

The estimates for x

and x

are more or less consistent but the estimate for x

does

differ, though all estimates of x

are very small.

The Brute force approach

Where possible, loops are to be avoided. In MATLAB this can be generally done

using vectorization. However, because of the “if statements”, this cannot be done

simply in the present case (it is possible to avoid looping by using ﬁlters but these

make very obscure code). Therefore, we shall employ a loop structure, which in

this circumstance runs fast enough not to warrant a more reﬁned approach (clear

code is to be preferred unless it impedes speed excessively). The program works in

a similar fashion as the R program in generating three vectors, x

, x

, and W and

then using a MATLAB function, max, to locate the maximum value of W and the

associated values of x

and x

function W¼FITNESS(v) % Function to evaluate fitness

{ This function is the same as that used in plotting}

Main Program:

clear all; % Clear the workspace

n¼100; % Number of divisions to be made

x1 ¼ linspace(0, 5, n); % x1 from 1 to 5 length n

x2 ¼ linspace(0, 5, n); % x2 from 1 to 5 length n

m ¼ n^2;

W ¼ zeros(m,1); x¼W;y¼W; % Preallocate x,y, W vectors

% Fitness values at each x1 x2 coordinate

row¼0; % Set index value to

% Iterate over all combinations of x1 and x2

for i ¼ 1:n

for j ¼ 1:n

row¼rowþ1; % Increment row

x(row)¼ x1(i); y(row)¼x2(j); % store x1 and x2

FISHERIAN OPTIMALITY MODELS 159

W(row)¼-FITNESS([x1(i),x2(j)]); % Convert to positive fit-

ness & store

end

[C,I] ¼ max(W); % C ¼ max W, I ¼ row in vectors

X1 ¼ x(I); X2 ¼ y(I); % Get values of x1 and x2

% Calculate x3

R ¼ 400; N ¼ 100; % Parameter values

x3 ¼ (R-N*(X1þX2))/N; % x3

[X1,X2,x3,C] % print out x1, x2, x3 and Wmax

OUTPUT:

ans ¼ 2.7778 1.2121 0.0101 2.3877

Using the brute force results to reﬁne the search (replace with the following

lines):

x1 ¼ linspace(2.7, 2.8, n); % x1 from 1 to 5 length n

x2 ¼ linspace(1.0, 1.3, n); % x2 from 1 to 5 length n

gives

ans ¼ 2.7697 1.2303 0 2.3880

Of course the results are the same as obtain in R.

2.18.31 Scenario 13: Plotting the ﬁtness function

The ﬁtness function follows that given in the R code. It is possible to replace the

loops doing the recursion with vectorized code using the function filter but the

code is so obscure that I would recommend it only if there is a real saving in time,

which in this case there certainly isn’t (for an example of using filter to replace

a loop see the online MathWorks support function (http://ww.mathworks.com/

support/tech-notes/1100/1109.html, p. 8).

function W¼FITNESS(x) % Function to evaluate fitness given Alpha

G ¼ x(1); % G values

Alpha ¼ x(2); % Alpha value

% Set parameter values

Mj ¼ 0.05; % Background mortality rate

Ma ¼ 0.4; % Constant of mortality function

Agemax ¼ 30; % Maximum age (arbitrary)

a ¼ 0.05; % Fecundity constant

A ¼ 10; % Wt increase/annum without reproduction

Aminus1 ¼ Alpha-1; % Year before first reproduction

S ¼ zeros(Agemax,1); % Pre-allocate vector for annual survival

% Growth prior to reproduction is linear

160 MODELING EVOLUTION

% To include cases in which growth is more complex I here use a for

loop

Wt ¼ zeros(Agemax,1); % Initialize Wt vector

S(1)¼ exp(-Mj); % Survival to age 1

% Calculate Wt and Survival from age 2 to alpha-1

for i ¼ 2: Aminus1

Wt(i) ¼ Wt(i-1)þ A; % Weight

S(i)¼ S(i-1)*exp(-Mj); % Annual survival

end

% Now calculate change in wt and survival for age alpha to max age

for i ¼ Alpha: Agemax

Wt(i)¼ Wt(i–1) þA-G*Wt(i–1); % Weight

S(i)¼ S(i–1)*exp(-(MjþMa*G)); % Annual survival

end

% Calculate W¼Fitness Note “.” for element by element vector multi-

plication

W ¼ a*sum(S(Alpha:Agemax).*Wt(Alpha:Agemax)*G);

W ¼W; % Return negative of fitness

Main Program:

clear all; % Clear the workspace

nG ¼ 11; % Number of G values

G ¼ linspace(0.01, 0.5, nG); % G from 0.01 to 0.5 length nG

nalpha ¼ 20-3þ1; % Number of ages to consider

alpha ¼ linspace(3, 20, nalpha); % x2 from 1 to 5 length 20

z ¼ zeros(nG,nalpha); % Preallocate z matrix

% Fitness values at each G alpha coordinate

for i ¼ 1:nG

for j ¼ 1:nalpha

z(i,j)¼-FITNESS([G(i),alpha(j)]); % Convert to positive fitness

end

subplot(2,2,1) % Divide graph sheet into 2 x 2 panels, contour

top left

[C,h]¼contour(alpha,G,z); % Create contour plot

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

lines.

clabel(C,h);

xlabel(’Alpha’); ylabel(’G’); % Add text

subplot(2,2,2) % Divide graph sheet into 2 x 2 panels, 3D top right

[xx,yy] ¼ meshgrid(alpha,G); % Create grids for 3D plot

surfc(xx,yy,z); % Plot 3D

xlabel(’Alpha’); ylabel(’x2’); zlabel(’Fitness’) % Add text

FISHERIAN OPTIMALITY MODELS 161

2.18.32 Scenario 13: Finding the maximum using a numerical approach

Brute force using many values

Because we want to calculate the optimal G for a given Alpha, we must make a

slight modiﬁcation to the ﬁtness function:

Replace

function W¼FITNESS(x) % Function to evaluate fitness given Alpha

G ¼ x(1); % G values

Alpha ¼ x(2); % Alpha value

with

function W¼FITNESS(G,Alpha) % Function to evaluate fitness given

Alpha

The above changes could have also been made in the plotting section. The

ﬁtness function is called by the main program:

clear all; % Clear the workspace

% Create vector for alpha values

nalpha ¼ 20-3þ1; % Number of ages to consider

alpha ¼linspace(3, 20, nalpha); % alpha vectorfrom 1 to 5 length 20

G ¼ zeros(nalpha,1); % Create vector to store best G values

W ¼ zeros(nalpha,1); % Create vector to store W values

for i ¼ 1:nalpha % Iterate over alpha vector values

% Find best G for this alpha by calling fminsearch

% alpha(i) is passed as a fixed parameter

G(i) ¼ fminsearch(@(G) FITNESS(G,alpha(i)),0.1); % Store best G

W(i) ¼ -FITNESS(G(i),alpha(i)); % Get W at best G for given alpha

end

% Now locate best combination and write out values

[C,I] ¼ max(W); % C ¼ max W, I ¼ row in vectors

Alpha_best ¼ alpha(I); % Best alpha

G_best ¼ G(I % Best G

W_best ¼ C; ); % W at best alpha, best G

[Alpha_best,G_best,W_best] % Print out best alpha, G and max W

OUTPUT: Ans ¼ 8.0000 0.1345 2.6849

Brute force using iteration

The ﬁtness function is the same as in the previous code. We also need another

function which I shall call BESTG that calculates the optimal G value for consecu-

tive pairs of Alpha. Note that this function passes back a vector.

% Function to get best G for consecutive pairs of alpha

function Wdiff ¼BESTG(alpha)

% Results for alpha

G1 ¼ fminsearch(@(G) FITNESS(G,alpha),0.1); % Store best G

W1 ¼ -FITNESS(G1,alpha); % Fitness

% Results for alphaþ1

162 MODELING EVOLUTION

G2 ¼ fminsearch(@(G) FITNESS(G,alphaþ1),0.1); % Store best G

W2 ¼ -FITNESS(G2,alphaþ1); % Fitness

W3 ¼ W2-W1; % Diff between fitnesses

% return Wdiff, W1, G1 % G1 will eventually be the best G

Wdiff ¼ [W3, W1, G1];

Main Program:

clear all; % Clear the workspace

ALPHA ¼ 5; % Set initial alpha

DIFF ¼ BESTG(ALPHA); % Calculate difference between W at two alphas

while DIFF(1)> 0; % If DIFF[1] > 0 then W still increasing

ALPHA ¼ ALPHAþ1; % Increment alpha

DIFF ¼ BESTG(ALPHA); % Call BESTG and get difference

end

% Out of loop and thus ALPHA is the best

[ALPHA, DIFF(3),DIFF(2)] % Print out alpha, G, W

2.18.33 Scenario 14: Finding the maximum using a numerical approach

MATLAB code changes in bold:

function W¼FITNESS(G,Alpha,N) % Function to evaluate fitness

given Alpha

% G IS NOW A VECTOR WITH N ENTRIES

% Set parameter values

Mj ¼ 0.05; % Background mortality rate

Ma ¼ 0.4; % Constant of mortality function

Age_max ¼ Alphaþ1; % Maximum age

a ¼ 0.05; % Fecundity constant

A ¼ 10; % Wt increase/annum without reproduction

A_minus_1 ¼ Alpha-1; % Year before first reproduction

S ¼ zeros(Age_max,1); % Vector of annual survival

% Growth prior to reproduction is linear

% To include cases in which growth is more complex I here use a for

loop

Wt ¼ zeros(Age_max,1); % Initialize Wt vector

S(1)¼ exp(-Mj); % Survival to age 1

% Calculate Wt and Survival from age 2 to alpha-1

for i ¼ 2 : A_minus_1

Wt(i) ¼ Wt(i-1)þ A; % Weight

S(i) ¼ S(i-1)*exp(-Mj); % Annual survival

end

% Now calculate change in wt and survival for age alpha to max age

% Accumulate W. W is now accumulated in the loop rather than after

W ¼ 0;

for J ¼ 1:N % Iterate over the adult ages

FISHERIAN OPTIMALITY MODELS 163

i ¼ Alpha þ J - 1; % Get age i ¼ Alpha, Alphaþ1

Wt(i)¼ Wt(i-1)þ A- G(J)*Wt(i-1); % Wt

S(i) ¼ S(i-1)*exp(-(MjþMa*G(J))); % Annual survival

W ¼ W þ a*S(i)*Wt(i)*G(J) ; % Cumulative fitness

end

W¼ -W; % Return negative of fitness

Now function BESTG:

% Function to get best G for consecutive pairs of alpha

function Wdiff ¼ BESTG(alpha)

N ¼ 2; % Nos of mature ages

% Results for alpha Note that G is passed two values

G1 ¼fminsearch(@(G) FITNESS(G,alpha,N),[0.1,0.1]); % Store best G

W1 ¼ -FITNESS(G1,alpha,N); % Fitness

% Results for alphaþ1

G2 ¼fminsearch(@(G)FITNESS( G,alphaþ1,N),[0.1,0.1]);% Store best G

W2 ¼ -FITNESS(G2,alphaþ1,N); % Fitness

W3 ¼ W2-W1; % Diff between fitnesses

% return Wdiff,W1,G1 % G1 will eventually be the best G

Wdiff ¼ [W3, W1, G1];

Main program calls BESTG:

clear all; % Clear the workspace

ALPHA ¼ 5; % Set initial alpha

DIFF ¼ BESTG(ALPHA); % Calculate difference between W at two alphas

while DIFF(1)> 0; % If DIFF[1] > 0 then W still increasing

ALPHA ¼ ALPHAþ1;

DIFF ¼ BESTG(ALPHA);

end

% Out of loop and thus ALPHA is the best

% Print Alpha, Wdif, G1, G2 ...GN

[ALPHA,DIFF(1),DIFF(3),DIFF(4)] %Print out alpha, G, W

OUTPUT: (Same as R)

Ans ¼ 19.0000 00.0003 0.3497 0.4835

To increase the number of ages and hence the number of variables we alter

N in

function BESTG,

and modify the following lines given that we wish to set N ¼ 5

(giving 6 variables to be estimated):

N ¼ 6; % Nos of mature ages

G1 ¼fminsearch(@(G) FITNESS(G,alpha,N), [0.1,0.1,0.1,0.1,0.1,0.1]);

G2 ¼fminsearch(@(G) FITNESS(G,alphaþ1,N),[0.1,0.1,0.1,0.1,0.1,0.1]);

[ALPHA,DIFF(1),DIFF(3),DIFF(4),DIFF(5),DIFF(6)] %Print out

The result is

ans ¼ 18.0000 0.0042 0.2166 0.2463 0.2917 0.3666

164 MODELING EVOLUTION

CHAPTER 3

Invasibility Analysis

3.1 Introduction

An alternative approach to that used in the last chapter is invasibility analysis,

which consists of asking if a clone displaying an alternate life history can

invade a resident population. While one could compare results for markedly

different life histories, in general, invasibility analysis has been used to locate

the optimal combinations of parameter values rather than qualitatively

different life histories. As with the “Fisherian” optimality approach, sexual repro-

duction is ignored. Invasibility analysis is used extensively, and is most

useful, when ﬁtness is density-dependent and there is age- or stage-structuring

in the model. The method can handle stable, cyclical, and chaotic population

dynamics. In this section I ﬁrst consider the general structure of age- and stage-

structured models and then describe the two general approaches of invasibility

analysis, namely pairwise-invasibility and multiple-invasibility analysis. For all

that you ever wanted to know about matrix population models see Caswell (2002).

3.1.1 Age- or stage-structured models

Consider the life table shown in Table 3.1.

Table 3.1 Calculation of age ‐speciﬁc survival probabilities and fertilities for the Leslie matrix.

Age xl

Post‐breeding census Pre‐breeding census

x−1

¼ S

xþ1

¼ S

010–– 0.8 0

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

1 0.80 1 0.80 0.8 0.4 0.4

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

2 0.20 3 0.40 1.2 0.25 0.75

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

3 0.05 4 0.25 1.0 0 0

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

4 0.00 1 0.00 0.00 ––

Using the post-breeding census, which is an assumption for the models discussed

in Chapter 2, the number of individuals entering age 1 at time t þ 1 is given by

1;tþ1

¼ S

1;t

þ S

2;t

þ S

3;t

þ S

4;t

¼ð0:8Þð1Þn

1;t

þð0:4Þð3Þn

2;t

þð0:25Þð4Þn

3;t

þð0Þn

4;t

¼ 0:8n

1;t

þ 1:2n

2;t

þ 1n

3;t

þ 0n

4;t

ð3:1Þ

where n

i,t

is the number in age class i at time t. The number surviving from t to

t þ 1 is given by

2;tþ1

¼ S

1;t

¼ 0:8n

1;t

3;tþ1

¼ S

2;t

¼ 0:4n

2;t

4;tþ1

¼ S

3;t

¼ 0:25n

3;t

ð3:2Þ

This set of equations can be represented in matrix format as

1;tþ1

2;tþ1

3;tþ1

4;tþ1

000

0 S

00S

1;t

2;t

3;t

4;t

ð3:3Þ

If a pre-breeding census is assumed, then S

¼ l

xþ1

(Caswell 1989, p. 12). The

matrix can be written in shorthand as

tþ1

¼ An

ð3:4Þ

The matrix A is known as the Leslie matrix after the ecologist who ﬁrst intro-

duced it into the population biology literature (Leslie 1945). The advantage of

using the matrix notation is that there are well-deﬁned rules for manipulating

matrices, particularly for matrix multiplication. From an evolutionary biologist’s

point of view the important feature of this matrix is that the population rate of

increase at a stable age distribution, l, is given by the ﬁrst eigenvalue of the

matrix. This value is readily obtained in R or MATLAB. For example, in the life

history speciﬁed by equations (3.1) and (3.2) the Leslie matrix is

A ¼

0:81:21:00

0:80 0 0

00:400

000:25 0

ð3:5Þ

which can be entered using R as

Leslie.matrix <- matrix(c(0.8, 1.2, 1.0, 0,

0.8,

0.0, 0.0, 0,

0.0, 0.4, 0.0, 0,

0.0, 0.0, 0.25, 0 ),4,4, byrow=TRUE)

where, for ease of writing, I have aligned the columns. The eigenvalues and

eigenvectors

can be obtained with the command eigen (R) or eig (MATLAB).

166 MODELING EVOLUTION

In R the appropriate eigenvalue can be drawn from the list with the sufﬁx

$values[1]. Thus the following commands in R,

Eigen.data <- eigen(Leslie.matrix)

Eigen.data$values[1] # Get first eigenvalue

gives 1.5516.

Equation (3.3) deﬁnes the change in population size after one generation. If the

initial population consists of a single gravid female the population size in the next

generation is given by

0:81:21:00

0:80 0 0

00:400

000:25 0

0:8

ð3:6Þ

Matrix multiplication in R is coded by %*%, thus

n <- Leslie.matrix%*%n

gives the multiplication shown in equation (3.6). Progressive application of

matrix multiplication produces the population projection. The following

coding calculates and plots the trajectories of the individual cohorts (ages 1–4)

and the total population size. Additionally, the observed rate of increase, given

by N

tþ1

, and the instantaneous rate of increase, r, given by log

tþ1

) are

plotted.

rm(list¼ls()) # Clear workspace

Leslie.matrix <- matrix(c(0.8, 1.2, 1.0, 0,

0.8, 0.0, 0.0, 0,

0.0, 0.4, 0.0, 0,

0.0, 0.0, 0.25,0),4,4, byrow¼TRUE)

Eigen.data <- eigen( Leslie.matrix)

Lambda <- Eigen.data$values[1] # Get first eigenvalue

Maxgen <- 12 # Number of generations

simulation runs

n <- c(1,0,0,0) # Initial population

# Pre-assign matrix to hold cohort number and total population size

Pop <- matrix(0,Maxgen,5)

Pop[1,] <- c(n[1:4], sum(n)) # Store initial population

# Pre-assign storage for observed lambda

Obs.lambda <- matrix(0,Maxgen,5)

for ( Igen in 2:Maxgen) # Iterate over generations

{

n <- Leslie.matrix%*%n # Apply matrix multiplication

INVASIBILITY ANALYSIS 167