Allman E.S., Rhodes J.A. Mathematical Models in Biology: An Introduction

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

76 Linear Models of Structured Populations

three-quarters of the adults die (after producing young) at each time

step, while the rest survive.

a. Model this situation using a matrix. Is this a Leslie or an Usher

model, or neither?

b. Compute eigenvectors and eigenvalues of the projection matrix

using MATLAB.

c. What is the intrinsic growth rate? The stable stage distribution?

2.3.10. Show that the absolute value for complex numbers satisﬁes

|(a + bi)(c + di)|=|a + bi||c + di|.

Projects

1. Consider a speciﬁc Leslie model with two age groups. After interpreting

each matrix entry in biological terms, investigate the behavior of your

model numerically using MATLAB for a variety of initial populations,

including the eigenvectors of the matrix. Explain how the eigenvalues

and eigenvectors are reﬂected in the behaviors you see when you plot the

populations over time. Repeat for several other matrices.

Suggestions

Begin with the Leslie model





using a MATLAB command sequence like:

P=[1/8 6; 1/5 0]

x=[10; 990]

xhistory=x

x=P*x, xhistory=[xhistory x]

plot(xhistory')

For a variety of choices of the initial populations, describe what ap-

pears to be happening to the populations over time. Do the number

of individuals in each group get bigger or smaller? Do they oscillate?

Compute the ratio of immature individuals to adults at various times.

How does this ratio change? Repeat this work with several different

2.3. Eigenvectors and Eigenvalues 77

choices of an initial vector. Qualitatively describe all the behaviors you

see.

Compute the eigenvectors and eigenvalues of A by entering [S,D]=

eig(A). Use the ﬁrst eigenvector as your initial vector by letting

x=S(:,1) and repeat the work above, including producing a plot.

Do this again using the second eigenvector with x=(:,2). Describe

the behavior of the model with these choices on initial vectors. How

is the behavior different? How is it the same? How are the eigenvalues

responsible for these behaviors?

How are the behaviors you see for the eigenvectors reﬂected in the

behavior you saw for other initial vectors?

Repeat all of the above for a few other models such as:





















Explain in biological terms why each of these models produces the

behavior it does. Then explain in terms of the eigenvalues and eigen-

vectors of the matrix why the behavior occurs.

Characterize the possible behavior of these 2 × 2 matrix models in

terms of the sign and size of the eigenvalues.

2. Leslie and Usher models can be used to help design intervention strategies

to help declining populations recover. A well-known example of this

was the study in (Crouse et al., 1987) on sea turtle populations that

supported the U.S. government-mandated use of turtle exclusion devices

in shrimpers’ nets.

An intervention might be designed to affect any one of the entries in the

Leslie matrix modeling a population. Because the dominant eigenvalue of

the matrix determines the overall growth rate, it is necessary to study how

changing each matrix entry affects the dominant eigenvalue. Determining

the effect of small changes in each of the entries is sometimes called

a sensitivity analysis. Imagine an endangered population grouped into

immature and mature subpopulations and modeled by the Usher model

with matrix



01.7

.3 .1



Analyze the effect of small changes in each of the nonzero entries in the

matrix on the population’s future.

78 Linear Models of Structured Populations

Suggestions

What is the dominant eigenvalue of the model? How fast will the

population increase or decline if no changes are made?

For the matrix



01.7

c .1



what values of c give a biologically meaningful model? For a variety of

values of c in that range, compute the dominant eigenvalue λ

. Present

the results of your computations in a table and as a graph of c vs. λ

Useful MATLAB commands are

lambda1vec=[]

cvec=[0:.1:1]

for c=cvec

A=[ 0 1.7;c .1]

lambda1=max(eig(A))

lambda1vec=[lambda1vec, lambda1]

end

plot(cvec, lambda1vec)

If you have read ahead to the next section of this text, ﬁnd a formula for

as a function of c? Does its graph agree with the one you produced?

If an intervention strategy attempts to change c in this matrix, describe

in biological terms what its focus might be. What value of c must be

achieved so that the population will recover?

Repeat the same sort of analysis to understand the affect of changing

the other nonzero entries of the matrix.

Without regard to the cost of implementing any recovery plan, which

entry do you think would be most effective to try to change? What bio-

logical issues might you need to understand better about the population

to answer this question adequately?

Why might a plan to change a fertility rate by a small amount have

different costs than a plan to change a survival rate?

Perform a sensitivity analysis on a Leslie or Usher model described by

a larger matrix of your choice.

2.4. Computing Eigenvectors and Eigenvalues

We ﬁrst show how eigenvectors and eigenvalues can be computed by hand

for 2 × 2 matrices.

2.4. Computing Eigenvectors and Eigenvalues 79

Given a matrix A, the eigenvector equation we wish to solve is Ax = λx,

where both the vector x and the scalar λ are unknown. This equation can be

rewritten as

Ax − λx = 0,

Ax − λI x = 0,

(A −λI)x = 0.

Notice that we had to stick the identity matrix into the middle equation

so that factoring the x out of each term would be valid. Without the identity,

we would have had A − λ, which makes no sense since we never deﬁned

subtraction of a scalar from a matrix.

Now, if x and λ are really an eigenvector and its eigenvalue, this last equa-

tion shows that the matrix (A − λI ) cannot have an inverse. For if it did,

then we could multiply each side of ( A − λI)x = 0 by the inverse to get

x = ( A − λI )

−1

0. But even without knowing what (A −λI)

−1

is explicitly,

we can tell that (A − λI )

−1

0 = 0, and this would mean x = 0. But our deﬁ-

nition of eigenvectors requires that they be nonzero, so we know this cannot

happen.

Now, if (A − λI ) does not have an inverse, then det(A − λI ) must be 0.

We have shown then that if λ is any eigenvalue of a matrix A, it must satisfy

the equation

det(A −λI) = 0.

Example. For the forest model matrix

P =



.9925 .0125

.0075 .9875



, P − λI =



.9925 − λ.0125

.0075 .9875 − λ



and so det(P − λI ) = 0 becomes

(.9925 − λ)(.9875 − λ) − (.0125)(.0075) = 0

− 1.98λ + .98 = 0

(λ − 1)(λ − .98) = 0.

This means the only possible eigenvalues for P are the numbers 1 and .98.

The equation det(A − λI ) = 0 is called the characteristic equation of A.

For a 2 × 2 matrix, it will always be a quadratic equation, so solving it will

be no harder than factoring or using the quadratic formula to ﬁnd the roots.

80 Linear Models of Structured Populations

While our argument above also applies to larger matrices (provided you

learn to compute determinants of larger matrices in some other course), solv-

ing the characteristic equation will be much harder because for an n × n

matrix, it involves an nth degree polynomial. Unless n is quite small, or you

are very lucky, this is simply not practical to do by hand. Nonetheless, we can

see that there can be at most n eigenvalues, since the characteristic equation

can have at most n roots. We also see that complex numbers may enter into

our calculations since the roots of a polynomial may not be real.

In summary, we have:

Theorem. If λ is an eigenvalue for an n × n matrix A, then it satisﬁes the

nth degree polynomial equation det( A − λI ) = 0. Thus, there are at most n

eigenvalues for A.

Once we have determined possible eigenvalues for a matrix, we need

to ﬁnd corresponding eigenvectors. We illustrate the process for P =



.9925 .0125

.0075 .9875



and λ = .98. We need to ﬁnd a vector v such that (P −

.98I )v = 0, so we must solve



.9925 − .98 .0125

.0075 .9875 − .98











.0125 .0125

.0075 .0075









Since we cannot solve this by ﬁnding the inverse of the matrix (why not?),

we write out the two equations this represents in nonmatrix form.

.0125x + .0125y = 0,

.0075x + .0075y = 0.

While guessing a nonzero solution is a perfectly valid way to proceed, we’ll be

more methodical in our approach. Because one of these equations is a multiple

of the other, to solve them we just need to solve .0125x + .0125y = 0. With

only one equation in two unknowns, we can take one of the unknowns to

have any value we like, and let that determine the second. For instance, if

we solve for y in terms of x,weﬁndy =−x. Thus, any vector of the form

v = (x, y) = (x, −x) = x(1, −1) is an eigenvector with eigenvalue .98. Since

we have freedom to choose x as we wish, we’ll take it to be 1. Thus, we have

found the eigenvector (1, −1) that was used throughout this chapter.

2.4. Computing Eigenvectors and Eigenvalues 81

The eigenvector associated with λ

= 1 can be found similarly:

P − 1I =



.9925 − 1 .0125

.0075 .9875 − 1





−.0075 .0125

.0075 −.0125



so we must solve

−.0075x + .0125y = 0

.0075x − .0125y = 0.

Because the equations are multiples of each other, we solve −.0075x +

.0125y = 0 to get y =

.0075

.0125

x =

x,sov =





= x





. Choosing

x = 5 (since it makes the vector have the simplest form), we ﬁnd v = (5, 3).

Although this was only one example of calculating an eigenvector for a

particular matrix P, for any 2 × 2 matrix the procedure works the same way.

Although we will not prove it here, you will always ﬁnd one of the equations

is a multiple of the other and then be able to solve for y in terms of x (or x in

terms of y) to ﬁnd all the eigenvectors.

As with eigenvalues, calculating eigenvectors for 3 × 3 or larger matrices

is done analogously to the 2 × 2 case, although some additional complications

come in. We’ll leave a discussion of those to a full course in linear algebra

and instead suggest that you let MATLAB do the computations for you.

Computer methods of calculation. Actually, MATLAB and other com-

puter packages do not really calculate eigenvectors and eigenvalues in the

way described previously. Because the computation of these is so important,

not only for biological models but for a host of problems throughout science

and engineering, quite clever and sophisticated methods have been developed

and incorporated into many standard software packages.

Although we will not really explain any methods these packages use, we

will give a hint at one type of approach by discussing the power method.

Given A, pick any initial vector x

and compute x

= Ax

. According

to the Strong Ergodic Theorem, if λ

is the dominant eigenvalue of A with

corresponding eigenvector v

, then we should expect

to be closer to v

than x

was. Because we do not yet know what λ

is, we have to somehow

adjust x

to account for the growth factor. One way of doing this is to simply

divide each entry of x

by its largest entry to get a new vector we call x



This means x



will have one entry that is a 1 and will be “closer” to being an

eigenvector than x

was.

We can then repeat the process using x



in place of x

to get an even better

approximate eigenvector. Of course, we should then repeat the process again,

82 Linear Models of Structured Populations

and again, until we ﬁnd our approximate eigenvectors are not changing by

much.

Example. Suppose A =





and we choose x





. Then















≈



.33

1.00







.33

1.00





3.33

7.33





7.33



3.33

7.33



≈



.45

1.00







.45

1.00





3.45

8.18





8.18



3.45

8.18



≈



.42

1.00







.42

1.00





3.42

7.96





7.96



3.42

7.96



≈



.43

1.00







.43

1.00





3.43

8.01





8.01



3.43

8.01



≈



.43

1.00



Thus, v

≈



.43

1.00



is the dominant eigenvector, at least to within a digit

or so of accuracy. To get the dominant eigenvalue λ

, we just note





.43

1.00





3.43

8.01



= 8.01



.43

1.00



so λ

≈ 8.01.



Compute the eigenvalues and eigenvectors for this matrix exactly and

see if they agree with this result.

2.4. Computing Eigenvectors and Eigenvalues 83

While this technique has only given us the dominant eigenvector and eigen-

value, variations on the idea can ﬁnd others.

Problems

2.4.1. For A =



.9 .3

.1 .7



, B =





, and C =



−13



a. Find the eigenvalues by solving the characteristic equations.

b. Find an eigenvector for each eigenvalue.

2.4.2. For the matrices in the preceeding problem, use the power method

to ﬁnd the dominant eigenvectors and eigenvalues. Do your answers

agree with those you found before?

2.4.3. Use the power method to ﬁnd the dominant eigenvector and eigenvalue

of the Leslie matrix





0073

.04 0 0

0 .39 .65





Check your answer by asking MATLAB for the eigenvectors and eigen-

values.

2.4.4. Actually, things can be more complicated than the text may have led

you to believe:

a. For A =





, ﬁnd two eigenvalues and two eigenvectors. (They

really do exist.)

b. For B =





, ﬁnd two eigenvalues and try to ﬁnd two eigenvec-

tors. What goes wrong?

2.4.5. Explain the connection between the power method and the Strong

Ergodic Theorem for linear models.

Nonlinear Models of Interactions

Our attention so far has been focused on modeling single populations. Al-

though we have broken a single population into subgroups, such as by age or

developmental stage, we have still treated it as if it is unaffected by the other

species or populations with which it might share an environment. Although

these models have provided valuable insights into how population sizes can

change, we now move our attentions to interactions between species or pop-

ulations.

Most living things interact with many coinhabitants of their environment.

Preying on other species, whether plant or animal, is a common way of taking

in energy; and most organisms are at risk of being preyed on themselves. But

not all important interactions between species are so obvious. Species may

ﬁnd themselves in competition for limited resources, whether food or space,

so that growth in one population is detrimental to another. Mutualism, where

several species interact in a way that beneﬁts both, also occurs in nature. A real

ecosystem may have hundreds or thousands of interacting populations, with all

sorts of direct and indirect interactions between them. How can we understand

the effects of these interactions without being lost in the complexity of it all?

To begin to understand the dynamics of interacting populations in such

systems, we start by focusing on only two populations and a single type of

interaction. The questions that will guide us are modiﬁcations of those we

have already asked when modeling a single population. For instance, what

mathematical formulas might model an interaction such as that of a predator

and its prey? What behavior does a computer simulation of such a model

show? Does one species disappear, and if so why? Do the populations reach

some equilibrium, oscillate, or jump wildly? Can such a system of interacting

populations show stability, and if so, under what circumstances?