• Stable age distribution and reproductive value
• Sensitivity analysis (optional for first go around)

Demography of the vegetable ivory palm Pytelephas

seemanii in Colombia, and the impact of seed harvesting

by Prof. Matt Miller, Department of Mathematics, University of South Carolina

email: miller@math.sc.edu

Maple worksheet to accompany Rodrigo Bernal,1998 ; Journal of Applied Ecology 35 :64-74.

This paper is an analysis of a stage-based model of the dynamics of a palm tree species, using, a modification of the Leslie model. In this case, the assumption is that the stage of an organism (size class) is a better indicator of reproductive performance and survival than age is. This assumption works well for many plants, and any organism in which size is correlated with survival and fecundity.

The general form of the model treats the stages as elements of a vector. The population projection matrix A (sometimes called a Lefkovich matrix, or, among mathematicians, transition matrix) is used to calculate the number of individuals in each stage class at time t+1, based on the numbers at time t.

n( t+1 ) = A n( t )

The (non-zero) elements of the matrix A are:

Pi the probabilities of surviving and remaining in the same stage.

Gi the probabilities of surviving and moving to the next stage.

Fi the fertility, or number of female offspring at time t+1, per adult female at time t.

> restart: with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> A:=matrix(4,4, [ [0,F2,F3,0], [G1,P2,0,0], [0,G2,P3,0], [0,0,G3,P4] ]);

For test purposes, try Brault & Caswell 1993 data. Confirm that all results are the same as published in that paper.

> Caswell:=matrix(4, 4,[ [0,0.0043,0.1132,0], [0.9775,0.9111,0,0],
[0,0.0736,0.9534,0], [0,0,0.0452,0.9804] ] );

Stable age distribution and reproductive value

The transition matrix is defined in Table 3 on p 70.

> A:=matrix(6,6,[[0.7052,0,0,14.8173,14.9770,14.7732],
[0.1548,0.6339,0,0,0,0],
[0,0.0216,0.9100,0,0,0],
[0,0,0.0330,0.9614,0,0],
[0,0,0,0.0386,0.9220,0],
[0,0,0,0,0.0188,0.9535]]);

>

The elements on the top row are the reproductive outputs. The subdiagonal elements represent the probability of transition to the next stage. The diagonal elements represent the probability of remaining in a stage.

***What fraction of seeds remain in the 'seed bank' each year?

***Of 1000 stage 1 individuals (seeds), how many survive to become stage 2? If there are currently 1000 stage 2 and 1000 stage 3, how many seeds (stage 1) will be present after 1 year? Suppose an individual is in the highest stage ). What is the probablility that this individual will still be alive 14 time steps from now?***

> eigsys:=eigenvects(A); # (output printing suppressed)

Out of the complete list of eigenvalues we find the maximum eigenvalue (lambda = population growth rate), and the index to the corresponding eigenvector (stable age distribution). This method for finding the population growth rate is described at the bottom of p 280..

> EVec:= seq( eigsys[i][3], i = 1 .. nops([eigsys])); # eigenvectors

> EVal:= seq( eigsys[i][1], i = 1 .. nops([eigsys])); # eigenvalues

> EM:= seq( eigsys[i][2], i = 1 .. nops([eigsys])); # multiplicities

> mags:=op( map( abs, [EVal] )); # magnitudes of eigenvalues

>

Next we identify the dominant eigenvalue (lambda), and the eigenvector it belongs to.

The dominant eigenvalue is called the finite rate of increase in the last paragraph of p 282.

> lambda:=max(mags);

> member( lambda, [mags], 'position'); indx:= position;

> V:=op( EVec[indx] );

The reported estimate of lambda in the paper is in Figure 3. The original transition matrix is for a population with no harvesting. How does the Maple estimate compare to the value in the Figure?

>

If the population grows for many generations with constant transition probabilities, the population stage vector will approach the eigenvector associated with the dominant eigenvector of the original matrix A.

This is the stable stage distribution vector, calculated from the dominant eigenvector of the original matrix A. Typically the stable stage vector w is scaled so that the sum of the elements is 1.0. To scale the stable stage vector, we have to do some fiddling. Add up the elements in the vector; then divide all elements in the vector by this sum:

> total:= sum( V['j'] , 'j' = 1 .. rowdim(A) );

>

We now have the stable stage distribution vector:

> w:=scalarmul(V, 1/total);

>

The stable stage distribution is reported as column SSD in Table 5 on p 70. How does it compare to our estimates?

>

Let's check; is A w = (lambda) w?

Remember this is our definition of an eigenvalue and eigenvector - that the eigenvector multiplied by the matrix is the same as the eigenvector multiplied by the eigenvalue. If we got the correct eigenvalues and eigenvectors, we should obtain the same answer by multiplying the stable age distribution by matrix A or by the eigenvalue:

> evalm( A &* w ); scalarmul( w , lambda);

Here we transpose the A matrix in order to set up calculation of the reproductive value vector, which is the dominant eigenvector of the transposed matrix.

> AT:=transpose(A);

> eigsysT:=eigenvects(AT): # (output printing suppressed)

Find maximum eigenvalue and index to corresponding eigenvector; this eigenvector is the reproductive value distribution, i.e., the left eigenvectors of matrix A ( last paragraph on p68)..

> EVecT:= seq( eigsysT[i][3], i = 1 .. nops([eigsysT]) ); # eigenvectors

> EValT:= seq( eigsysT[i][1], i = 1 .. nops([eigsysT]) ); # eigenvalues

> EMT:= seq( eigsysT[i][2], i = 1 ..nops([eigsysT]) ); # multiplicities

> magsT:= op( map( abs , [EValT] )); # magnitudes of eigenvalues (Note: [ ] turns a sequence into a list; op turns a list into a sequence.)

> mu:= max( magsT); # get dominant eigenvalue

> member(mu, [magsT], 'position'); indexT:= position;

Convert eigenvector lists to actual vectors for use in sensitivity analysis.

> VT:= op( EVecT[indexT] );

This is the reproductive value distribution. The reproductive values are reported as columns of Table 5 on p 70. Note that vector v is scaled so the first element is 1.0. We do a scalar multiplication of all elements in the vector by the reciprocal of the first element. (Recall that a scalar multiple of an eigenvector is again an eigenvector with the same eigenvalue.) We now have the reproductive value vector corresponding the to the column "v", printed in Table 5 on p 70:

> v:= scalarmul(VT, 1/VT[1] );

>

How do these values compare to those printed in the table? Remember once again that the authors of the paper probably used a different numerical accuracy in their calculations.

Let's confirm that we correctly identified the eigenvalue/eigenvector pair.

Let's just check that AT v = (mu) v.

> evalm( AT &* v ); scalarmul( v , mu);

>

Sensitivity analysis (optional for first go around)

Now we will calculate the sensitivity of the system to changes in fecundity. This is the analysis proposed in Caswell. We are analytically determining the partial derivative of lambda (population growth rate) with respect to each element in the projection matrix. Thus we are asking: given a change in a particular element in the projection matrix, what is the corresponding change in population growth rate? This allows us to determine which components of the life history have the biggest effects on population growth, and which may be most sensitive to natural selection. The first parameter we need is the inner product of the two eigenvectors (stable age distribution and reproductive value distribution). This is symbolized as <w,v> in the Caswell handout. This value is used as the denominator of the calculation.

> IP:= innerprod(VT , V); # dot, or scalar, product

> sensitivity_to_fecundity:= [seq( V[i]*VT[1] / IP , i = 1 .. rowdim(A) )];

>

Observe that it is immaterial whether we use the normalized vectors or not, as long as we do this consistently.

> [seq( w[i]*v[1] /innerprod( w, v) , i = 1 .. rowdim(A) )];

In terms of the paper this corresponds to the top row of matrix S, equation 7, in Caswell except that the first and last elements have been suppressed because the original matrix had zero entries in these locations.

We do a similar calcuation to obtain the sensitivity to survival:

> sensitivity_to_survival:= [seq( V[i]*VT[i+1] / IP , i = 1 .. rowdim(A)-1 )];

This vector corresponds to the subdiagonal of matrix S. The entire sensitivity matrix can be calcuated by generalizing the above operations. Note there are some entries that can be taken to be zero, or otherwise ignored, because the A matrix had zero entries (but this isn't really necessary, since as we shall see in a moment we end up multiplying each entry of S by the corresponding entry of A).

> S:=matrix(rowdim(A),rowdim(A), [ seq( seq( V[i] * VT[j] / IP, i = 1.. rowdim(A)), j = 1.. rowdim(A) ) ]);

The elasticity values are calculated by scaling the sensitivities by aij / lambda (Caswell handout). The resulting matrix is reported as Table 4 on p 70. The elasticities are a bit easier to interpret than the sensitivities, because they represent the proportional contributions of matrix elements to lambda. Recall the we can't use the label E in Maple because this means the base of natural logarithms (2.718...).

> e:=matrix(rowdim(A),rowdim(A),[ seq( seq( S[j,i] * A[j,i] / lambda, i = 1 .. rowdim(A)), j = 1.. rowdim(A)\)]);

***Try seeing if you believe that the sensitivities or elasticities are correct, following the pattern illustrated below.*** We will change one of the entries in the original "A" matrix by 10% and see what percent change there is in lambda. Compare this change with the corresponding value in the elasticity or sensitivity matrix. Let's first recall the original matrix, eigenvalues, eigenvectors:

> A:= evalm(A);

> lambda[old]:= lambda; mu[old]:= mu;

> RV:= evalm(v); SAD:= evalm(w);

Let's change the P2 value by10% (we'll do this by multiplying P2 by 1.1) and see what percent change we get in lambda. This is the so-called [2,2] entry of the matrix A, that is, the number that appears in the first row (from top to bottom) and second column (from left to right). If the author's description of elasticity is correct, then we should get a change in lambda that is proportional to the elasticity value for P2.

> A[2,2]:=1.1*A[2,2]; # make the change in the affected entry

> A:=evalm(A); # confirm the matrix is OK

Now let's calculate the eigenvalues and eigenvectors for this revised matrix.

> eigsys:=eigenvects(A): # (output printing suppressed)

> EVec:= seq( eigsys[i][3], i = 1 .. nops([eigsys])); # eigenvectors

> EVal:= seq( eigsys[i][1], i = 1 .. nops([eigsys])); # eigenvalues

> mags:= op( map( abs, [EVal] )); lambda[new]:=max(mags);

Compare this new value for lambda with the old value of lambda:

> lambda[old];

This is the ratio of the new value to the old; the fractional change in lambda is the portion after the decimal. Did you expect lthe new lambda to be larger? Why?

> ratio:= lambda[new]/lambda[old];

Here is the elasticity value. Compare it to the fractional change in lambda.

> elas:= e[2,2];

Note that we made a 10 percent change in one element of the matrix (a fractional change of 0.1). If you multiply this elasticity value by the fractional change made in the matrix element, you should get a number close to the fractional change in lambda. If you made a 100% change (a fractional change of 1.0) in the matrix element, you should see a fractional change in lambda equal to the elasticity value.

> (elas * 0.1) - (ratio - 1); # close to zero?

> member( lambda[new], [mags], 'position'); indx:= position;

> evalm(SAD); # here's the original stable age distribution

Compare with the new stable age distribution.

> V := op( EVec[indx] ); total:= sum( V['k'] ,' k' = 1 .. rowdim(A) );
w:= scalarmul(V, 1/total); # compare with SAD above

Transpose the A matrix in order to set up calculation of the reproductive value vector, which is the dominant eigenvector of the transposed matrix.

> AT:=transpose(A);

> eigsysT:=eigenvects(AT): # (output printing suppressed)

> EVecT:= seq( eigsysT[i][3], i = 1 .. rowdim(A) ); # eigenvectors

> EValT:= seq( eigsysT[i][1], i = 1 .. rowdim(A) ); # eigenvalues

> magsT:= op( map( abs , [EValT] )); # magnitudes of eigenvalues

> mu[new]:= max( magsT); member(mu[new], [magsT], 'position'); indexT:= position;

> VT:= op( EVecT[indexT] ); v:= scalarmul(VT, 1/VT[1] ); # compare with RV above

> evalm(RV);

>

>

>

***Now, how did v and w change, given that 10% change in a single entry of the matrix A? Considering the change that was made can you give an explanation for the result? Now change the survival probability G3 by 15% (multiply it by 1.5). (Don't forget to reset A[1,2] back to its original value first) What effect do you expect to see? What happens in fact?***

>

The author states on p71 that "an 86% reduction in fecundity of all classes is required to reduce lambda from 1.0589 to 1.0". Try making reductions in fecundity and see if your observations confirm this: First we will get an original copy of matrix A

> A:=matrix(6,6,[[0.7052,0,0,14.8173,14.9770,14.7732],
[0.1548,0.6339,0,0,0,0],
[0,0.0216,0.9100,0,0,0],
[0,0,0.0330,0.9614,0,0],
[0,0,0,0.0386,0.9220,0],
[0,0,0,0,0.0188,0.9535]]);A := matrix([[.7052, 0, 0, 14.8173, 14.9770, 14.7732], [.1548, .6339, 0, 0, 0, 0], [0, .216e-1, .9100, 0, 0, 0], [0, 0, .330e-1, .9614, 0, 0], [0, 0, 0, .386e-1, .9220, 0], [0, 0, 0, 0, .188e-1, .9535]]);

>

Then we define our fractional reduction in fecundity.

> fecundity_reduction:=0.86;

We will now reduce the seed production entries by this amount. We will not change A[1,1], because that element of the matrix defines the rate of retention of seeds in the seed bank, not a fecundity. Therefore we change only elements 2 to 6 in the top row of the matrix:

> for i from 2 to rowdim(A) do

> A[1,i]:=A[1,i]*(1-fecundity_reduction); od:

> evalm(A);

Now move back to the top of the worksheet, and rerun from a point below the original definition of matrix A. See if the values for lambda, stable stage distribution, and reproductive value distribution change as expected.

>