Some simple models of cannibalism

Some Simple Models of Cannibalism JAMES

C. FRAUENTHAL

Department

of Mathematics,

Harvey Mudd College, Claremont,

Crrlifornia 9171 I

Recemed I June 1982: revised 4 August 1982

ABSTRACT A sequence of mathematical models for a species which engages in cannibalism are investigated. The models treat the species as age-structured, and assume that adults consume the unhatched eggs of their own kind. The McKendrick or von Foerster partial differential differential dynamical

equation model is first converted into a set of three coupled, nonlinear ordinary equations, and then adjusted to describe cannibalism. Some rather unusual effects are discovered. These include both a Hopf bifurcation and a catastrophic

transition

from an asymptotically

stable equilibrium

point to a stable limit cycle.

INTRODUCTION Several papers have appeared recently which investigate predator-prey interactions in which the prey species is age-structured and the predators consume only the eggs of the prey. This paper extends this idea by replacing the predator species with the adults of the prey species. In other words, it investigates mathematical models of cannibalism in which the adults of a species consume the eggs of their own kind. The first paper to appear on egg-eating predators is by Gurtin and Levine [l] and investigates a model analogous to the well-known Lotka-Volterra equations. Since the effect of egg predation appears to be a destabilizing influence, it is not surprising that the ordinarily neutrally stable LotkaVolterra equations become unstable. The authors report growing amplitude oscillations which lead ultimately to the extinction first of the prey species and then (since no alternative food source is explicitly included in the model) of the predator species. In a sequel to this paper, Levine [2] introduces several more-realistic effects into the model. In particular, the author finds that introducing either a carrying capacity for the prey or mutual interference between the predators leads to a stable equilibrium with both species present. The equilibrium in the model with the prey carrying capacity can be either an asymptotically stable point or a limit cycle, while the model with predator interference appears to lead just to an asymptotically stable point. MATHEMATICAL

BIOSCIENCES

63:87-98

OElsevier Science Publishing Co., Inc., 1983 52 Vanderbilt Ave., New York, NY 10017

87

(1983)

0025.5564/83/01087+

12$03.00

88

JAMES

C. FRAUENTHAL

Thompson, Di Biasio, and Mendes [3] investigate the same model that Levine [2] calls predator interference, though they arrive at the formulation by a very different route than Levine. These authors also conclude that the system has an asymptotically stable equilibrium point with both species present. Coleman and Frauenthal[4] study a model with both a prey carrying capacity and predator appetite satiation. Their equilibrium point with both species present experiences a Hopf bifurcation at a critical parameter value and manifests peculiarly shaped limit cycles. Most of the previous models for cannibalism have been rather complicated in structure, and have appeared in the resource management literature. The interested reader should see, for example, McKelvey, Hankin, Yanosko, and Snygg [5] and Botsford and Wickham [6], two papers concerned with recently observed oscillations in the number of adult Northern California Dungeness crabs (Cancer magister), a species known to practice cannibalism. Although the authors disagree as to the cause of the cycles, they do both reproduce a figure based upon catch data indicating large amplitude, periodic oscillations of the crab species. A recent paper by Gurtin and Levine [7] discusses models for cannibalism which are similar to the ones in the present paper. Both papers consider phenomenological models, without attempting to fit parameters to data. The present paper considers two points not looked at by Gurtin and Levine. Specifically, it investigates the stability of periodic solutions around a unique positive equilibrium, and considers the case of multiple positive equilibria. In the process, a mechanism is described which could account for the sudden onset of cyclic behavior.

THE MATHEMATICAL

MODEL

We begin with the well-known McKendrick [8] or von Foerster equation for the evolution of an age-structured population:

[9]

where p( a, t) is the age density (the number of individuals of age a at time t), B(t) is the egg production rate, P(t) is the total population, ~(a, P) is the mortality function, and P(a, P) is the maternity function. As shown by

SOME SIMPLE

MODELS

Gurtin and MacCamy

89

OF CANNIBALISM

[lo], if it is assumed that /J(a, p) = P(P) /3(a,P)

only, (2)

lY> 0,

=P(P)aepa”,

then the partial differential equation can be converted, by finding moments with respect to a weighting function g(u), to a closed set of three nonlinear ordinary differential

equations:

P=-/.l((P)P+p(O,t)

]g=ll,

B=-[p(P)+cu]B+p(P)A

[g=aeKau],

k=-[/.l(P)+“]A+p(O,t)

[g=

where ’ denotes time differentiation the relation

A(f) =/,,e

(3)

eda’],

and A is an auxiliary variable defined by

-““p(a,t)

da.

We take (3) as the starting point for our model. Note that in arriving at these equations we have assumed that the mortality is age-independent (not a very good assumption) and that the maternity as a function of age just rises and then falls (a realistic assumption for many species). We proceed now to specialize the equations further. In doing so we must come to grips with the distinction between B(t) and ~(0, r). Since B(t) represents the egg production rate and ~(0, t) the number of newborns, their difference is the number of unhatched eggs. We will assume that all unhatched eggs have been the victims of cannibalism by adults of the species being modeled. Since it is impossible for the adults to survive exclusively by cannibalism, we assume that there are alternative food sources available which are not explicitly included in the model. Hence, the number of eggs which are eaten depends upon the product of the number of eggs and the number of adults, but with a satiation factor included to account for the fact that in times of large egg production, the adults do not modify their customary diet in direct proportion with egg availability. Such a mechanism has been discussed by Holling [ 111, though the original analytic form which will be used appears to be due to Michaelis and Menten [ 121. Based upon this idea, we represent the number of cannibalized eggs as

B(t)-p(O,t)=

Nt)P(t) 1-t kB(t)

Strictly, there should be one more multiplicative

’

k>O

(4)

constant

on the right hand

90

JAMES C. FRAUENTHAL

side; however, without loss of generality it can be absorbed into the parametrization. The consequence of this reparametrization is to scale P, B, and A all by the same constant. Thus, while the relative sizes of the three dependent variables remain significant, their absolute sizes are unimportant. Solving (4) for ~(0, t) immediately presents a small difficulty. It is formally possible for the number of newborns to become negative; consequently we write p (0, t) = max

B -

(5)

To complete the model we assume that mortality is a linearly increasing function of population size (an assumption which leads to a carrying capacity for the species), but that the birth rate is population-independent, so

Combining

(3), (S), and (6) leads to the set of equations:

(7)

where all of the parameters are positive. Before continuing it is worth noting that the maximum function is inconvenient to treat analytically. Throughout this paper, when equations are solved analytically, the nonzero term in the maximum function is used, and then checked a posteriori for suitability. In numerical simulations, the maximum function is employed, since it presents no computational difficulties.

THE MODEL WITHOUT CARRYING OR APPETITE SATIATION

CAPACITY

Although we have set out to caricature reality with a simple model, the set of equations we have derived do not yield easily to analysis. We therefore proceed initially to suppress certain of the features we have included in the model by setting some parameters to zero. Later we will put the suppressed features back into the formulation, but as we do so we will be forced to move more to numerical methods.

SOME SIMPLE

MODELS

91

OF CANNIBALISM

We begin by suppressing both the carrying capacity (EL,= 0) and the appetite satiation of the adults (k = 0), and consider the model P=-j~aP+max{B-BP,O}, B=-(&+++&A, k = -(pO This set of equations librium point at

(p*,B*,A*)=

(8)

+ OL)A+max{B

- BP,O}.

is easily shown to have its only nontrivial

[ I-

equi-

(~o~a~2]( 1,(plfo)2 &i).

We restrict our attention to the situation PO > (p. + LU)~,as this corresponds to positive numbers of animals. So long as this inequality holds, it is easily shown that the trivial equilibrium point at the origin is unstable. If, however, the inequality is reversed, the equilibrium at the origin (extinction) becomes stable. It is interesting to note that the analogous model with cannibalism omitted leads to unbounded exponential growth; hence cannibalism introduces an equilibrium point. Linearized stability analysis of the nontrivial equilibrium point leads to a cubic characteristic equation all of whose roots can be shown using the Routh array to lie in the left half of the complex plane. This means that the equilibrium point is stable with respect to small perturbations. Using both numerical root extraction on the characteristic equation and integration of the differential equations (8) it was concluded that the solutions always approach the equilibrium (and do so in an oscillatory manner only when the mortality parameter, po, assumes very small values).

THE MODEL WITHOUT

CARRYING

CAPACITY

We next extend the model by reintroducing the appetite satiation term; hence we make k > 0. This means that the governing equations now assume the form P=-p,P+max

A=-(p,+a)A+max This set of equations

i

B-&,0

(

>

B-&,0

,

I

is easily seen to have just one nontrivial

equilibrium

92

JAMES C. FRAUENTHAL

point. It occurs at

(p*,B*,A*)=

l-(Po+42/Po 1+kpo[l-&/(p0+a)2]

, i

POP0 ‘tPo+ny~oP~”

1’

which is just a constant scaling of the equilibrium point found for (8). Once again, to ensure that the equilibrium point corresponds to positive numbers of animals, we require p > (p. + a) 2. Linearized stability analysis proceeds exactly as in the earlier case except that it is no longer possible to show that all three roots of the characteristic equation lie in the left half of the complex plane. Instead of attempting to continue analytically, we employ three separate numerical methods. These require that we select numerical values for the four parameters. These have been chosen so as to illustrate the existence of a Hopf bifurcation. We choose p. = 0.015, cr = 1.0, k = 1.5, and continue to treat PO as a parameter. The first numerical procedure simply evaluates the roots of the characteristic equation associated with the linearized stability analysis of the nontrivial equilibrium point. Numerical experimentation confirms that when (p. + a) 2 < Do < 2.83, the one real and two complex conjugate roots have negative real parts, whereas when PO > 2.84 (though less than about 8.3), the complex conjugate roots have positive real parts, though the real root remains negative. This suggests that somewhere on the interval 2.83 < PO < 2.84, the local properties of the equilibrium point switch from stable to unstable as PO is increased. In order to confirm whether a limit cycle will develop, we next employ a numerical routine which tests for the existence of a Hopf bifurcation. The software package used is fully described in a book by Hassard, Kazarinoff, and Wan [ 131, and was graciously provided by the first author. Applying the routine to our equations confirms that indeed a Hopf bifurcation into a stable limit cycle occurs at PO = 2.8368. As a final test, we solved the set of differential equations (9) using a very fast and efficient solver [ 141 which is available on the Harvey Mudd College computer. Computer simulations with parameter choices on either side of the predicted critical parameter for the Hopf bifurcation are shown in Figures 1 and 2. With reference to the Northern California Dungeness crab population, this appears to represent an incorrect mechanism to explain the onset of oscillatory behavior following exploitation. Since this species is harvested by removing adult males, one might expect the effect of the fishery to be a reduction in the value of the fertility as measured by PO, since fewer males would make it harder for females to find mates.

SOME SIMPLE

0

MODELS

100

200 model

FIG. I. capacity

Computer

93

OF CANNIBALISM

simulation

300

400

500

time

of the model with appetite

saturation

but no carrying

(9) with p0 = 0.015, a= 1.0, k = 1.5, and p,, = 2.0.

THE FULL

MODEL

We proceed finally to study the full model as stated in (7). Before doing so it is worth noting that if we include a carrying capacity (CL,> 0) but no appetite saturation (k = 0), we arrive at results which look almost exactly like those for the first case considered. The existence of the carrying capacity adds a second nontrivial equilibrium point, but it always corresponds to negative numbers of animals. We therefore do not dwell on this case, but move directly to the full model. The set of equations for the location of the nontrivial equilibrium points is now a quartic. Due to the difficulty of doing general analysis on such an

JAMES

n

C. FRAUENTHAL

r

r B

0

100

200 model

: : : , ! J: 300

; : :L: 400

:/, 500

time

FIG. 2. Computer simulation of the model with appetite capacity (9) with p0 = 0.15, a = 1.0, k = 1.5, and PO = 2.84.

saturation

but no carrying

equation with five parameters, we immediately fix three of the parameters at the same values as in the previous example (pe = 0.015, (Y= 1.0, and k = 1.5) and set )J, = 0.07. We look at results as a function of the value of &,. It appears that for sensible values of &, there are always at least two real equilibrium points, one with positive and one with negative numbers of animals. As & is decreased from about 5.0 to about 4.8, two new real, positive equilibrium points are born; one then migrates toward the previously existing, real, positive equilibrium point; as these two points coalesce, they annihilate, ceasing to exist as real roots. This situation is shown in Figure 3, which can be visualized as a “time lapse” picture taken as &, is varied.

SOME SIMPLE

MODELS

OF CANNIBALISM

,\

-1 0

4

2 population

6 size

, 8

P

FIG. 3.

capacity varying

The population size equilibrium points for the model with both carrying and appetite saturation (7) with p,, = 0.015, p, = 0.07. a = 1.0, k = 1.5, and PO

from 4.8 to 5.0. The symbols

S and U denote,

respectively,

that the equilibrium

point is locally stable or unstable.

What are particularly interesting are the local stability properties of the three real, positive equilibrium points. Numerical stability analysis as in the previous cases confirms that the largest equilibrium point is always locally asymptotically stable and nonoscillatory, while the other two equilibrium points are both locally unstable. The middle one is nonoscillatory, and the smallest one is oscillatory. The phenomenon that we observe as PO decreases from about 5.0 to about 4.8 is often called a catastrophe. Imagine for example that & = 4.9. The only stable, positive equilibrium point is nonoscillatory, so the population would be expected to remain at this equilibrium and manifest no tendency to cycle.

96

JAMES

C. FRAUENTHAL

A computer-generated trace of this situation is shown in Figure 4. Then suppose that some external effect shifts &, downward slightly to PO = 4.8. The stable nonoscillatory equilibrium point vanishes, and the system shifts catastrophically to the only remaining positive equilibrium position. But around this equilibrium point there already exists a strongly attractive limit cycle. This solution is illustrated in Figure 5. Although catastrophic shifts between locally stable equilibrium points are not particularly rare, the shift from a stable point to a stable limit cycle is unusual. This mechanism might be closely related to the one which caused the population of the Northern California Dungeness crab to suddenly start oscillating. As noted at the end

B

6

L_ P

, 0

:

., 3 00

:

: : / : : : : I 300 200 model

:

:

:

:

/

:

400

:

:

:

I

500

time

FIG. 4. Computer simulation of the model with both carrying capacity saturation (7) with p0 = 0.015, w, = 0.07, (Y= 1.0, k = 1.5, and & = 4.9.

and appetite

SOME SIMPLE

MODELS

97

OF CANNIBALISM

P

\

\ \

1 \

300

200

mode

I i-

400

FIG. 5. Computer simulation of the model with both carrying capacity saturation (7) with p,, = 0.015. p, = 0.07, OL= 1.0, k = 1.5, and &, = 4.8.

of the previous section, in the Northern

500

time

California

Dungeness

and appetite

crab fishery,

adult males are harvested, likely causing a reduction in the value of &,. The full model manifests a catastrophic cycle as & is reduced.

CONCLUDING

shift from a stable point to a stable limit

REMARKS

A mathematical model has been proposed to explain some of the phenomena observed for species which practice cannibalism. The model is analyzed using mostly numerical techniques. For this reason, the analysis is frequently

98

JAMES

C. FRAUENTHAL

specialized to regions of parameter space in which interesting (though not necessarily common or realistic) phenomena occur. In summary, we have found the following results. If the model excludes both a carrying capacity and adult appetite satiation, the solution is a single, locally stable equilibrium point. Adding just a carrying capacity has little effect, while adding just appetite saturation permits a single equilibrium point to experience a Hopf bifurcation into a stable limit cycle. However, this bifurcation is associated with increasing fertility, a situation opposite to the one which may have led to the onset of cycles in the Northern California Dungeness crab fishery. Finally, adding both a carrying capacity and appetite satiation allows the model to admit a catastrophic change from a locally stable equilibrium point to a fully developed limit cycle as fertility is reduced slightly. REFERENCES 1 2

M. E. Gurtin and D. S. Levine, On predator-prey interactions with predation dependent on age of prey, Math. Biosci. 47:207-219 (1979). D. S. Levine, On the stability of a predator-prey system with egg-eating predators,

3

Math Biosci. 56127-46 (1981). R. W. Thompson, D. Di Biasio, and C. Mendes, predators, Mafh Biosci., 60: 109- 120 (1982).

4

C. S. Coleman 63:99-l

5

6

and

J. C. Frauenthal,

Satiable

Predator-prey egg-eating

interactions: predators,

Egg-eating Math.

Biosci.

19 (1983).

R. McKelvey. D. Hankin, K. Yanosko, and C. Snygg, recruitment models: An application to the Northern

Stable cycles in multi-stage California Dungeness crab

(Cancer magisfer) fishery, Canad. J. Fish. Aquut. Sri. 37:2323-2345 (1980). L. W. Botsford and D. E. Wickham, Behavior of age specific, density dependent models and the Northern California Dungeness crab (Cancer magister) fishery, J. Fish. Res. Bd. Canada, 35:833-843 (1978).

7

M. E. Gurtin and D. S. Levine, On populations Appl. Math. 42:94-108 (1982).

8

A. C. McKendrick, Applications Math. Sot. 44:98-130 (1926).

9

H. von Foerster, Some remarks on changing populations, in The Kinetics of CeNulcrr Proliferafion, Grune and Stratton, New York, 1959. M. E. Gurtin and R. C. MacCamy, Some simple models for nonlinear age-dependent population dynamics, Math. Biosci. 43 : 199-2 11 ( 1979).

IO II 12 13 14

of mathematics

that cannibalize to medical

their young,

problems,

SIAM

J.

Proc. Edinburgh

C. S. Holling, The functional response of predators to prey density and its role in mimicry and population regulation, Mem. Ent. Sot. Canada 45:3-60 (1965). L. Michaelis and M. L. Menten, Biochem. Z. 49:333 (1913). B. D. Hassard, N. D. Kazarinoff, and Y.-H. Wan, Theory and Application of Hopf Bifurcation, London Math. Sot. Lecture Notes X41, Cambridge U.P., 1981. N. Freed and K. Carosso, Omnigraph project report, Harvey Mudd College Software Library Guide, 198 I.