Age-structured Population Models—The Leslie Matrix



This worksheet requires some knowledge of Linear Algebra. The material needed for this worksheet is covered in Sections 9.2 and 9.3 of Neuhauser (2004) Calculus for Biology and Medicine. 2nd edition. Prentice Hall. A summary is provided below.

Matrix Algebra

A [pic] matrix [pic], is a rectangular array of numbers with m rows and n columns. We write this as

[pic]

A [pic] matrix is called row vector, and a [pic] matrix is called a column vector. Two [pic] matrices A and B are the same if all their entries agree. (We can only compare matrices of the same size.)

Matrices of the same size can be added. Suppose that A and B are [pic] matrices. Then

[pic]

is a [pic] matrix with entries

[pic]

for [pic] and [pic]. Matrix addition satisfies the following two properties:

[pic]

Matrices in MATLAB

To generate a vector in MATLAB, type, for instance,

>> u=[2 4 6 8 10]

This same vector can also be created by typing

>> u=2:2:10

If you want to square each entry, type

>>v=u.^2

Note the dot before the “^”. It indicates that each component is squared.

To enter a matrix, use semicolons. Entries within the same row may be separated by commas:

>>A=[1, 0, 3; 2, 4, 5]

To add two matrices, first define two matrices that are of the same size, for instance

>>A=[1,2,3;4,5,6]

>>B=[-1,2,-3;-4,5,-6]

and then add them

>>C=A+B

The transpose of a matrix A, denoted by A’, is obtained by switching rows and columns. If you type

>>D=A’

you can see the effect.

Matrix Multiplication

Matrix multiplication is quite a bit more involved. We give the definition first. Suppose [pic] is a [pic] matrix and [pic] is a [pic] matrix, then

[pic]

is a [pic] matrix with

[pic]

for [pic] and [pic]. Note that in this product, the number of columns of the matrix A must be the same as the number of rows of the matrix B.

Here is an example. Assume that

[pic] and [pic]

Since the number of columns of A is the same as the number of rows of B, the product [pic] is defined, and C is a [pic] matrix. To find the entry in say the second row and the third column of the matrix C, you would multiply the numbers of the second row of matrix A by the numbers of the third column of the matrix B.

Matrix Multiplication in MATLAB

Enter the matrices A and B and the multiply them:

>>A=[1,2,3;-1,0,4]

>>B=[1,2,3,-3;0,-1,4,0;-1,0,-2,1]

>>C=A*B

You will see that c23 is indeed -11.

To get a better understanding of how MATLAB deals with multiplication of matrices, enter the following matrix

>>D=[1 2;3 4]

This is a [pic] matrix, so the product DD is defined.

>>D*D

This is the same as

>>D^2

Make sure you distinguish this from the following two operations that use a dot before the operator. The first one multiplies each of the corresponding entries (and is thus different from matrix multiplication):

>>D.*D

Likewise,

>>D.^2

takes the square of each entry but does not square the matrix.

Age-structured Population Models—The Leslie Matrix

In many species, reproduction is highly age-dependent. For instance, periodical cicadas spend 13-17 years in the nymphal stage; they only reproduce once in their lifetime. Many animals, such as humans, elephants, etc., do not reproduce during their first years and then their reproductive success is age-dependent. To model such situations, age-dependent population models are appropriate. Patrick Leslie introduced matrix models that have discrete age classes with synchronous reproduction (Leslie 1945). The models are parameterized by age-specific survival probabilities and average number of female offspring.

An age-structured population with k age classes is described by a population vector of length k

[pic]

where [pic] denotes the number of females in age class j. Individuals within an age class are assumed to have equal birth and death probabilities. Survival from age class j to j+1 is given by probability Pj; female individuals in age class j give birth on average to Fj female offspring. This can be illustrated using life cycle graphs (see Figure 1):

[pic]

Figure 1: A life cycle graph with four age classes

Translating the life cycle graph into equations, we find

[pic]

This can be written in matrix notation

[pic]

The matrix contains both the survival probabilities (just below the diagonal) and the fecundities (in the first row). The matrix is called Leslie matrix. Take a moment to verify that the two descriptions are equivalent.

Task 1:

Suppose a Leslie matrix is of the form

[pic]

(a) Translate the matrix into a life cycle graph. (b) If at time 0, there are 100 individuals in age class 1, and no individuals in each of age classes 2 and 3, how many individuals are in the three age classes at time 2? You may use MATLAB to solve this. Note that if you enter the population vector at time 0 as a row vector, you need to transpose that vector first before you can compute the population vector at time 1.

To analyze matrix models, we need to learn some tools from linear algebra.

Eigenvalues and Eigenvectors

Definition: Assume that A is a square matrix. A nonzero vector w that satisfies the equation

[pic]

is a (right) eigenvector of the matrix A and the number λ is the corresponding eigenvalue.

Note that we require that the eigenvector is different from the zero vector. An eigenvalue can be 0, it can even be a complex number. The geometric interpretation of eigenvectors explains its importance. An eigenvector corresponding to a given matrix is a vector that gets mapped by the matrix on the line through the eigenvector. The eigenvalue allows us to find the length of the mapped vector.

To find eigenvectors and eigenvalues, we need to solve the equation

[pic]

We can rewrite this as

[pic]

where I is the identity matrix. (The identity matrix is a square matrix that has 1s on the diagonal line and 0s elsewhere.) In order to obtain a nontrivial solution (i.e., [pic]), the matrix [pic] must be singular. The determinant of a matrix can be used to check whether a matrix is singular, and the condition is:

[pic]

Solving this equation will yield values for λ. For each value of λ, we can then find an eigenvector w. Determinants were originally introduced to determine whether a system of linear equations has a unique solution. The theory of determinants has moved far beyond this. But to us, this will be the situation where we will encounter determinants.

If A is a two-by-two matrix, the computation can be carried out by hand. Namely, if

[pic]

then

[pic]

Multiplying this out yields a quadratic equation, which can be solved using

[pic]

If all eigenvalues of the matrix A are distinct and real, the eigenvectors will be real and can be computed by solving [pic] for each eigenvalue.

Determinants in MATLAB

The determinant of a matrix can be computed in MATLAB. Define a square matrix A. Then the command is

>>det(A)

To use MATLAB to find the eigenvalues of a matrix using the determinant, the following commands are used (below, we’ll introduce a quicker way):

>>syms lam

>>A=[2 3;1 4];

>>I=[1 0;0 1];

>>solve(det(A-lam*I));

Eigenvalues and Eigenvectors in MATLAB

Suppose the matrix A is given as

[pic]

To compute the eigenvalues of A in MATLAB, there is a command that is much faster than the series of commands we used previously:

>>eig(A)

The answer is [pic] and [pic]. A calculation by hand confirms this, namely,

[pic]

Factoring the right-hand side and then setting it equal to zero yields

[pic]

and hence the same answer.

The eigenvectors can be computed by hand once the eigenvalues are known by solving the equation

[pic]

In the case of the eigenvalue [pic], this amounts to solving

[pic]

which results in the following system of linear equations

[pic]

This reduces to just one equation since the two equations are identical:

[pic]

This has infinitely many solutions. We need to pick one, for instance, the vector [pic] satisfies the equation.

To find the eigenvalues and eigenvectors using MATLAB, type

>> [V,D] = eig(A)

This produces two matrices. The matrix D is a diagonal matrix with the eigenvalues along the main diagonal. The columns of the matrix V are the eigenvectors scaled to length 1. In this case, we find

[pic] and [pic]

Note that the length of a vector [pic] is calculated as [pic].

Stable Age Distribution

In Step 2b of the Computer Lab, we will compute the fraction of individuals in each of the age classes. If these fractions stabilize over time, they represent the stable age distribution. A stable age distribution does not always exist. Suppose that [pic] is a stable age distribution, then applying the Leslie matrix to this vector, that is computing [pic], results in a new vector that can be scaled so that it is equal to [pic]. This means that [pic] is a right eigenvector. The scaling can be expressed as [pic], implying that [pic] is the corresponding eigenvalues.

We will now use an important result from linear algebra, which we cannot prove here. Assume that the eigenvalues of the Leslie matrix are real and distinct. Because of the assumption of distinct eigenvalues, the corresponding right eigenvectors are also distinct, and we can write any initial population vector n(0) as a linear combination of the right eigenvectors of the Leslie matrix

[pic]

where the coefficients ci are determined by the initial vector n(0). Applying the Leslie matrix to this, we see what the advantage of writing the initial vector in this way is

[pic]

If we compute the population size at time 2, we find

[pic]

and so on. We thus find

[pic]

We find that the population vector at time t is a weighted average of exponentials [pic] where the weights are determined by the initial population vector.

If all the eigenvalues are distinct, we can order them so that [pic]. We then find that for t large, the largest eigenvalue dominates in the sense that

[pic]

which we can also express as

[pic]

The eigenvector w1 corresponding to the largest eigenvalue [pic] is also the stable age distribution. The largest eigenvalues [pic] tells us how quickly the population grows for large times. We see that when [pic], the population will grow exponentially without bound, and when [pic], it will decline exponentially fast.

Some Generalizations

So far, we assumed that all eigenvectors are distinct. We wish to relate this property to biological properties of the population. We begin with a number of definitions.

A matrix is nonnegative if all its entries are nonnegative. A matrix is positive if all its entries are positive. Note that all Leslie matrices are nonnegative.

Irreducibility

A nonnegative matrix or a life cycle graph is irreducible if and only if every stage can be reached from every other stage. In terms of the life cycle graph this means that there is a path from every node to every other node.

Life cycle graphs with postreproductive age classes do not result in an irreducible matrix.

Primitivity

An irreducible nonnegative matrix A is primitive if there exists a [pic] so that [pic] is positive.

A loop is a path that has the same start and end point and passes through each node at most once.

A life cycle graph is primitive if it is irreducible and the greatest common divisor of the lengths of its loops is 1 (see Figure 2).

[pic]

Figure 2: The life cycle graph on the left is not primitive since there is only one loop and it is of length 3. The life cycle graph on the right is primitive since there are two loops, one of length 2, the other of length 3; consequently, the greatest common divisor is 1.

An irreducible nonnegative matrix A is primitive if its associated life cycle graph is primitive.

Case 1

If the Leslie matrix is irreducible and primitive, then

[pic] with [pic]

In this case,

[pic]

and the right eigenvector w1 corresponding to the largest eigenvalue is a stable age distribution (see Computer Lab, Steps 1-2).

Case 2

If the Leslie matrix has postreproductive classes, it is reducible. The Leslie matrix can be reduced into four submatrices

[pic]

The matrix B1 alone gives the growth rate and is unaffected by the post reproductive age classes (see Computer Lab, Step 3).

Case 3

If the Leslie matrix is irreducible and nonprimitive, then

[pic]

This corresponds to cyclic behavior (see Computer Lab, Step 4).

In-class Lab (due __________________)

As preparation for this lab,

• Read the worksheet carefully and do the tasks.

• Read the steps below.

• Read the article by R.N. Holdway and C. Jacomb. 2000. Rapid Extinction of the Moas (Aves: Dinornithiformes): Models, Test, and Implications.

Step 1

In the accompanying EXCEL spreadsheet (Tab: “Leslie Matrix”), the following Leslie Matrix is entered

[pic]

We will run the population for 50 time steps when initially there are 200 individuals in the first age class and no other individuals. Note that row 9 is different from the other rows (the other rows are explained at the end of this worksheet). In row 9, we entered the following in the respective cells:

Cell B9: =G3

Cell C9: =G4

Cell D9: =G5

Cell E9: =G6

Cell F9: =SUM(B9:E9)

Cell G9: =F10/F9

Cell H9: =B9/F9

Cell I9: =C9/F9

Cell J9: =D9/F9

Cell K9: =E9/F9

Row 10 is explained at the end of this worksheet. You can drag the rows down so that you can investigate the behavior for the first 50 time steps.

Step 2

(a) Plot the population sizes of each age class as a function of time for t=0,1,2,…,20. Log transform the vertical axis and replot. Explain why for sufficiently large t the data points lie along a straight line. Plot the proportion of each age class as a function of time and use it to find a stable age distribution. (Note that once you transform the vertical axis to a logarithmic scale, you might get a message from Microsoft Excel telling you that “Negative or zero values cannot be plotted correctly on log charts….” Simply change the scale back to the arithmetic scale and you get rid of these messages.)

(b) Column G lists the ratio of consecutive population sizes. What is the significance of this value? How is it related to the largest eigenvalues?

(c) Columns H-K list the fraction of individuals in each age class. Plot the fraction for each age class as a function of time for t=0,1,2,…,20. Is there a stable age distribution and if yes, what is it?

Step 3

Change the matrix from Step 1 so that there are two post reproductive age classes:

[pic]

Check that the growth parameter of the population is unaffected by changes in the survival probabilities of the post reproductive age classes by increasing P3 from 0 to 1 in increments of 0.1. Assume that initially there are 200 individuals in the first age class and no other individuals.

Step 4

Analyze the following nonprimitive Leslie matrix

[pic]

Assume that initially there are 200 individuals in the first age class and no other individuals. Describe what you see. Change the survival probabilities from 0.7 to 0.5. What happens?

Step 5

Steps 1-4 is done more easily in Matlab. Here is the code for Steps 1 and 2:

clear

L=[1.6 1.5 0.25 0.1; 0.8 0 0 0; 0 0.5 0 0; 0 0 0.3 0];

X0row=[200 0 0 0];

X0=X0row';

Time=101;

NPop=zeros(4,Time);

NPop(:,1)=X0;

for i=1:Time-1

NPop(:,i+1)=L*NPop(:,i);

end

Gen=0:Time-1;

semilogy(Gen,NPop); xlabel('Generation'); ylabel('Population Size');title('Leslie Model');

legend('Age 1','Age 2','Age 3','Age 4')

[V,D]=eig(L);

norm=sum(V(:,1));

age=V(:,1)./norm;

Modify the code so that you can do Steps 3 and 4 as well.

Step 6

Comment on the appropriateness of use of Leslie matrices in the article by Holdway and Jacomb in light of the fact that the Leslie matrix corresponds to exponential growth/decay.

Literature Cited

P.H. Leslie. 1945. On the Use of Matrices in Certain Population Mathematics. Biometrika

33:183-212.

  |A |B |C |D |E |F |G |H |I |J |K | |1 |Leslie Matrix | | | | | | | | | | |2 | |1 |2 |3 |4 | |Initial vector | | | | |3 |1 |1.6 |1.5 |0.25 |0.1 |N1= |200 | | | | | |4 |2 |0.8 |0 |0 |0 |N2= |0 | | | | | |5 |3 |0 |0.5 |0 |0 |N3= |0 | | | | | |6 |4 |0 |0 |0.3 |0 |N4= |0 | | | | | |7 | | | | | | | | | | | | |8 |Year |Stage 1 |Stage 2 |Stage 3 |Stage 4 |Total | λ_t |w_11 |w_12 |w_13 |w_14 | |9 |0 |200 |0 |0 |0 |200 |2.4 |1 |0 |0 |0 | |10 |1 |320 |160 |0 |0 |480 |2.266667 |0.666667 |0.333333 |0 |0 | |11 |2 |752 |256 |80 |0 |1088 |2.169853 |0.691176 |0.235294 |0.073529 |0 | |12 |3 |1607.2 |601.6 |128 |24 |2360.8 |2.174382 |0.680786 |0.254829 |0.054219 |0.010166 | |13 |4 |3508.32 |1285.76 |300.8 |38.4 |5133.28 |2.174198 |0.683446 |0.250475 |0.058598 |0.007481 | |

Cell B10: =$B$3*B9+$C$3*C9+$D$3*D9+$E$3*E9

Cell C10: =$B$4*B9+$C$4*C9+$D$4*D9+$E$4*E9

Cell D10: =$B$5*B9+$C$5*C9+$D$5*D9+$E$5*E9

Cell E10: =$B$6*B9+$C$6*C9+$D$6*D9+$E$6*E9

Cell F10: =SUM(B10:E10)

Cell G10: =F11/F10

Cell H10: =B10/F10

Cell I10: =C10/F10

Cell J10: =D10/F10

Cell K10: =E10/F10

-----------------------

F4

F1

F3

F2

P3

P2

P1

3

4

2

1

[pic]

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

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download