Pod-specific demography of killer-whales
by Prof. Matt Miller, Department of Mathematics, University of South Carolina email: miller@math.sc.edu
Maple worksheet to accompany Brault and Caswell, 1993
Pod-specific demography of killer-whales , Ecology 74 : 1444-1454.
This paper is an analysis of a stage-based model of the dynamics of killer whales, a modification of the Leslie model. In this case, the assumption is that the stage of an organism (yearling, juvenile, mature female, post-reproductive female) 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. This is equation 1 on p. 1445:
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] ]);
Stable age distribution and reproductive value
We will now define the matrix with the values listed on p. 1447:
> A:=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] ] );
This is equation 2 on p. 1447. 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. ***Of 1000 yearlings, how many survive to become juveniles? If there are currently 1000 juveniles (stage 2) and 1000 reproductive adults (stage 3), how many yearlings (stage 1) will be present after 1 year? Suppose an individual is in the highest stage (non-reproductive adult). 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 distribution is vector w in equation 5 on p. 1448.
> EVec:= seq( eigsys[i][3], i = 1 .. 4); # eigenvectors
> EVal:= seq( eigsys[i][1], i = 1 .. 4); # eigenvalues
> EM:= seq( eigsys[i][2], i = 1 .. 4); # multiplicities
> mags:=op( map( abs, [EVal] )); # magnitudes of eigenvalues
Next we identify the dominant eigenvalue (lambda), and the eigenvector it belongs to.
> lambda:=max(mags);
> member( lambda, [mags], 'position'); location:= position;
> V:=op( EVec[location] );
This is the stable age distribution vector, calculated from the dominant eigenvector of the original matrix A. Note that in equation 5, vector w is scaled so that the sum of the elements is 1.0. To scale the stable age 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 .. 4 );
We now have the stable age distribution vector as printed on p. 1448:
> w:=scalarmul(V, 1/total);
Let's check; is A w = (lambda) w?
> evalm( A &* w ); scalarmul( w , lambda);
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 vector "v" in equation 5 on p. 1448.
> EVecT:= seq( eigsysT[i][3], i = 1 .. 4 ); # eigenvectors
> EValT:= seq( eigsysT[i][1], i = 1 .. 4 ); # eigenvalues
> EMT:= seq( eigsysT[i][2], i = 1 ..4 ); # 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. Note that in equation 5, 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 as printed on p. 1448:
> v:= scalarmul(VT, 1/VT[1] );
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 on p. 1447 in equation 3. 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 equation 3 on p. 1447. This value is used as the denominator on the right hand size of equation 3. Observe that it is immaterial whether we use the normalized vecotrs or not.
> IP:= innerprod(VT , V); # dot, or scalar, product
> sensitivity_to_fecundity:= [seq( V[i]*VT[1] / IP , i = 1 .. 4 )];
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 .. 4 )];
In terms of the paper this corresponds to the top row of matrix S on p. 1448, equation 7, 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 .. 3 )];
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(4,4, [ seq( seq( V[i] * VT[j] / IP, i = 1.. 4), j = 1.. 4 ) ]);
The elasticity values are calculated by scaling the sensitivities by aij / lambda (equation 4 on p. 1447). The resulting matrix is reported as equation 8 on p. 1448. 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 E in Maple because this means the base of natural logarithms (2.718...).
> e:=matrix(4,4,[ seq( seq( S[j,i] * A[j,i] / lambda, i = 1 .. 4), j = 1.. 4)]);
***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 F2 value from 0.0043 to 0.00473 (a 10% change) and see what percent change we get in lambda. This is the so-called [1,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 Caswell's description of elasticity is correct, then we should get a change in lambda that is proportional to the elasticity value for F2.
> A[1,2]:=.00473; # 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.
> mags:= op( map( abs, [EVal] )); lambda[new]:=max(mags);
Compare this new value for lambda with the "estimate" of lambda at bottom of p 1447. This is 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[1,2];
Note that we made a 10 percent change in fecundity (a fractional change of 0.1). If you multiply this elasticity value by the fractional change made in the fecundity, 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 fecundity value, 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'); location:= position;
> evalm(SAD); # here's the original stable age distribution
Compare with the new stable age distribution.
> V := op( EVec[location] ); total:= sum( V['k'] ,' k' = 1 .. 4 ); w:= scalarmul(V, 1/total); # compare with SAD above
> 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% from 0.0452 to 0.05198. (Don't forget to reset A[1,2] back to its original value of 0.0043!) What effect do you expect to see? What happens in fact?***
Individual pod data
If you want to do a similar analysis on the individual pod data in the appendix, you need to set up a symbolic description of matrix A, and then substitute in the values from the table:
> Ag:=[[0,F2,F3,0],[G1,P2,0,0],[0,G2,P3,0],[0,0,G3,P4]];
> A:=matrix(4,4,subs({G1=0.9535,G2=0.0803,G3=0.0414,P2=0.8827,P3=0.9586,P4=0.9752,F2=0.0067,F3=0.1632},Ag));
Now move your cursor back to the line that does the first eigenvalue calculation at the top of the worksheet. This is the line that defines "eigsys" just after A is first defined. Then go back through the worksheet to get new eigenvalues, eigenvectors, sensitivities and elasticities. Alternatively, replace the lines defining A at the top of the worksheet with these lines, and use what has been a closely guarded secret up til now. In the menu bar on top look for Execute Worksheet (on some systems it's under View, but this varies). BAM!! DONE!!
>
NOTE - if you use new data in the worksheet, you may find that some of the procedures don't work. Usually the problems arise from specific matrices having repeated eigenvalues. We can show you how to work around this problem in class.
