Nonlinear Differential Equations, Electron. J. Diff. Eqns., Conf. 05, 2000, pp. 255–283 http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu or ejde.math.unt.edu (login: ftp)

On the complex formation approach in modeling predator prey relations, mating, and sexual

disease transmission ^{∗}

Horst R. Thieme & Jinling Yang

Dedicated to Alan Lazer on his 60th birthday

Abstract

Complex formation is used as a unified approach to derive representations and approximations of the functional response in predator prey relations, mating, and sexual disease transmission. Applications are given to the impact of a generalist predator on a prey population and the spread of a sexually transmitted disease in a multi-group heterosexual population.

1 Introduction

In mathematical models of such diverse bioscientific areas as demographics, ecology and epidemics the following quite different looking questions have complex (or pair) formation as a hidden common theme:

How does the number of marriages depend on the number of eligible male and female singles?

How does the number of prey which is killed and eaten depend on the number of prey available and on the number of predators?

How does the number of sexual contacts depend on the number of available individuals?

The correct modeling of pair formation (often called the two-sex problem) is an important issue in human demography, we refer to Hadeler et. al. (1988) for a discussion and the history of various marriage functions. Predator-prey relations are important in many ecosystems. The classical experimental and theoretical work by Holling (1965, 1966) studies the functional response of one predator to the amount of available prey, but does not necessarily give the functional response for several predators. The controversy concerning the paradox of enrichment has led to various modifications (Beddington, 1975, DeAngelis et al., 1975, e.g.), among others to so-called ratio-dependent models (see Arditi and Ginzburg (1989), Huisman and De Boer, 1997, Kuang and Beretta, 1998, Cosner et al., 1999, e.g., and the references mentioned there).

∗Mathematics Subject Classifications: 92D25, 92D30, 92D40.

Key words: Complex (pair) formation, functional response, approximations of Beddington,

Michaelis-Menten alias Monod alias Holling 2 and of Ross type, proportionate mixing,

generalist predator, disease incidence, basic replacement ratio.

\begin{matrix}c \\ ◯ \end{matrix}2000 Southwest Texas State University.

Published October 25, 2000.

H.R.T. was partially supported by NSF grants DMS-9403884 and DMS-9706787.

The question of the correct functional relationships becomes even more difficult, if several species are involved or if the population needs to be structured into various groups or according to age (see Castillo-Chavez and Velasco-Hernandez, 1994, for a survey and many references).

Once complex formation is revealed as underlying theme of these different phenomena, it is suggestive to apply the Michaelis-Menten (1913) quasi-steady-state approach from enzyme kinetics to model this process from first principles. Hsu and Fredrickson (1975) seem to be first to use this approach for marriage and mating models, Heesterbeek and Metz (1993) add mathematical rigor and include sexually transmitted diseases, and De Boer and Perelson (1995), Borghans, De Boer and Segel (1996) and Huisman and De Boer (1997) (and some earlier papers cited therein) apply it to T cell proliferation or/and predator-prey models.

Unfortunately an explicit solution of the quasi-steady state equations has been elusive in the multiple group or multiple species case. Hsu and Fredrickson (1975) find a closed solution in the case of homogeneous mixing, and De Boer and Perelson (1995), Borghans, De Boer and Segel (1996) and Huisman and De Boer (1997) present approximations in the case of one prey and many predator species or analogous situations.

The aims of this note are the following:

• present a unified modeling approach to mating functions, predator-prey relationships and contact rates in sexually transmitted diseases (Section 2 and Appendix A)

• derive approximations and estimates for the case of many prey and many predator species and analogous situations, including approximations of symmetric or asymmetric Beddington type, Michaelis-Menten alias Monod alias Holling II type, and of Ross type (Section 2 and Appendix B)

• use the explicit solution of the steady state equations in the case of proportionate mixing to compare the various approximations (Section 3)

• illustrate that interesting features of qualitative behavior may be lost in the impact of a generalist predator on a prey population, if the explicit solution of the quasi-steady-state approximation is replaced by approximations (Section 4)

• Use the complex formation approach to model the incidence function for sexual disease transmission in a multigroup hetero-sexual population, and illustrate how the basic reproduction number can be determined for the case of separable mixing in spite of seemingly overwhelming complexity (Section 5).

2 The formation of complexes

Assume that there are m species of prey and n species of predators and that X_{i} denotes the size of the i^{th} prey species and Y_{l} the size of the l^{th} predator species.

Alternatively let there be m groups of females and n groups of males in the sexually active part of a population with X_{i} denoting the sizes of the female groups and Y_{l} the sizes of the male groups. In this section these sizes are assumed to be constant in time.

Both predation and sexual contact require the formation of a complex which, for simplicity, consists of two individuals. We assume that the complex has two stages, engagement (courtship) and consumption (handling, mating).

Let C_{il} be the number of complexes in the engagement stage with one individual from X_{i} and the other individual from Y_{l}, while D_{il} is the number of complexes in the consumption stage. In a predator-prey model, D_{il} can also be interpreted as the number of predators of species l that are handling (eating and digesting) prey of species i (cf. Huisman, De Boer, 1997).

Let π_{j} be the average probability of an individual from X_{j} to be active, i.e., to be involved in a complex or to be available to form a complex, while π̃_{l} is the average probability of an individual from Y_{l} to be active. These numbers between 0 and 1 are called the activity or availability levels of the respective group or species. Since every active (or available) individual is either an active (or available) single or involved in a complex, we have the following laws,

π_{i}X_{i} = x_{i} + ∑_{k=1}^{n}[C_{ik} + κD_{ik}], i = 1, …, m,

(2.1)

π̃_{l}Y_{l} = y_{l} + ∑_{j=1}^{m}[C_{jl} + D_{jl}], l = 1, …, n,

where x_{i} is the number of active singles from X_{i} and y_{l} is the number of active singles from Y_{l}.

For a mating model κ = 1, while κ = 0 for a predator prey model because the prey is already dead in the consumption stage.

The dynamics of complexes is described by the following differential equations based on mass action type kinetics,

C′_{il} = γ_{il}x_{i}y_{l} − ρ_{il}C_{il},
(2.2)
D′_{il} = p_{il}ρ_{il}C_{il} − σ_{il}D_{il}.

Here 1/γ_{il} is the average length it takes a given pair of an available X_{i} single and an available Y_{l} single to get together, while 1/ρ_{il} is the average length of the first stage and 1/σ_{il} the average length of the second stage of the complex. Of course it is possible that a complex is also dissolved by natural death of one of the two individuals, but the duration of the stages is assumed to be so short that this can be neglected. p_{il} is the probability with which the complex actually progresses from the first to the second stage. In a predator prey model, this is the probability with which a predator from species i manages to kill a prey from species l without being killed itself in the process.

The following theorem concerning global existence and uniqueness of solutions of (2.1) and (2.2) is proved in Appendix A.

Theorem 2.1 For all initial data there exists a unique local solution of (2.1), (2.2). Solutions which are non-negative initially (meaning that also x_{i}, y_{l} are non-negative initially), exist for all forward times and remain non-negative.

At this point it must be mentioned that, in the applications we have in mind, equations (2.1) and (2.2) are only part of larger models and that X_{i} and Y_{l} vary in time, while we assume in this section that they are constant. One would like to simplify the large models by using a quasi-steady-state approach for (2.1), (2.2), which consists in replacing C_{il} and D_{il} by the steady states C_{il}^{∗} and D_{il}^{∗} of (2.1), (2.2). The basic assumption, which guarantees that one obtains a useful approximation this way, requires X_{i} and Y_{l} to vary much more slowly than the number of complexes. A rigorous mathematical justification can be based on the theory of singular perturbations. Heesterbeek and Metz (1993) prove for one stage complexes that one obtains useful approximations on finite time intervals. Conditions which guarantee the preservation of asymptotic stability can be derived from Hoppensteadt (1974).

In any case it is necessary that (2.1) and (2.2), for time-independent X_{i} and Y_{l}, have unique and globally asymptotically stable steady states. For one stage complexes this has been shown by Heesterbeek and Metz (1993) using Lyapunov functions. Using different methods, we at least show uniqueness (Appendix B) for the two stage model and give sufficient conditions for local and global asymptotic stability (Theorem 2.3, Appendices C and D). These conditions basically say that the time it takes to find a partner for forming a complex is long compared to the duration of the stages of the complex.

The equations for the steady state are obtained by setting the derivatives equal to 0 in (2.2),

C_{il}^{∗} = δ_{il}x_{i}^{∗}y_{l}^{∗}, D_{il}^{∗} = \frac{p_{il}ρ_{il}}{σ_{il}} C_{il}^{∗},

x_{i}^{∗} = π_{i}X_{i} − ∑_{k=1}^{n}[C_{ik}^{∗} + κD_{ik}^{∗}], (2.3)

y_{l}^{∗} = π̃_{l}Y_{l} − ∑_{j=1}^{m}[C_{jl}^{∗} + D_{jl}^{∗}],

with

δ_{il} = \frac{γ_{il}}{ρ_{il}} .  (2.4)

We expedite the algebra by introducing

x_{i}^{∗} = u_{i}π_{i}X_{i},
y_{l}^{∗} = w_{l}π̃_{l}Y_{l}, (2.5)
C_{il}^{∗} = K_{il}δ_{il}π_{i}X_{i}π̃_{l}Y_{l}.

Then

K_{il} = u_{i}w_{l},

u_{i} = 1 − ∑_{k=1}^{n} π̃_{k}Y_{k}ξ_{ik}δ_{ik}K_{ik}, (2.6)

w_{l} = 1 − ∑_{j=1}^{m} π_{j}X_{j}ξ̃_{jl}δ_{jl}K_{jl},

with

ξ_{jk} = 1 + κ\frac{p_{jk}ρ_{jk}}{σ_{jk}} , κ ∈ \{0, 1\},  ξ̃_{jk} = 1 + \frac{p_{jk}ρ_{jk}}{σ_{jk}} .  (2.7)

From (2.6) we can derive separate fixed point equations for the vectors u and w. In general it is not possible to solve these fixed point equations explicitly, but it can be shown that a unique solution exists and can be approximated by successive iterations (see Appendix B).

Theorem 2.2 There exists a unique equilibrium of the system (2.1), (2.2), i.e., a unique solution of the algebraic system (2.3). It can be found as limits of successive approximations for the vectors u and w in (2.6).

As Hsu and Fredrickson (1975) already noticed in an age-dependent mating model, an explicit formula for the equilibrium can be found under the assumption of proportionate mixing (see Section 3). In the general case, we at least have estimates from below and from above (Appendix B). In the following we write C and D rather than C^{∗} and D^{∗}. The estimate from below is reminiscent of Ross solutions (Castillo-Chavez, Busenberg, 1991),

C_{il} ≥ \frac{π_{i}X_{i}}{1 + ∑_{j=1}^{m} ξ̃_{jl}δ_{jl}π_{j}X_{j}} δ_{il} \frac{π̃_{l}Y_{l}}{1 + ∑_{k=1}^{n} ξ_{ik}δ_{ik}π̃_{k}Y_{k}} ,  (2.8)

the estimate from above has the following form,

C_{il} ≤ \frac{π_{i}X_{i}}{1 + ∑_{j=1}^{m} \frac{ξ̃_{jl}δ_{jl}π_{j}X_{j}}{1 + ∑_{k=1}^{n} ξ_{jk}δ_{jk}π̃_{k}Y_{k}}} δ_{il} \frac{π̃_{l}Y_{l}}{1 + ∑_{k=1}^{n} \frac{ξ_{ik}δ_{ik}π̃_{k}Y_{k}}{1 + ∑_{j=1}^{m} ξ̃_{jk}δ_{jk}π_{j}X_{j}}} . (2.9)

These estimates can be improved by iteration (Appendix B), but, as seen in the last one, reach the limits of manageability very fast.

If the γ_{il} are very small compared with the ρ_{il} and σ_{il} (e.g., if the time needed to find one prey is large compared with the time needed to eventually hunt it down, kill it, eat and digest it), formula (2.8) provides a reasonable approximation, which we call the Ross approximation, but we have reason to believe in general (and have proved it for proportionate mixing in Section 3) that extensions of the formulas derived by Beddington (1975) and DeAngelis et al. (1975) are as good or even better (Appendix B),

C_{il} ≈ \frac{δ_{il}π_{i}π̃_{l}X_{i}Y_{l}}{1 + ∑_{k=1}^{n} ξ_{ik}δ_{ik}π̃_{k}Y_{k} + ∑_{j=1}^{m} ξ̃_{jl}δ_{jl}π_{j}X_{j}} .  (2.10)

We mention that (2.10) is a lower estimate of C_{il} if m = 1 or n = 1 (Appendix B).

In prey-predator models where ξ_{ik} < ξ̃_{ik} with the difference being possibly considerable, it is suggestive to mix the estimates (2.8) and (2.9) in order to obtain an approximation (see Appendix B, in particular (B7) for the justification),

C_{il} ≈ \frac{π_{i}X_{i}}{1 + ∑_{j=1}^{m} ξ̃_{jl}δ_{jl}π_{j}X_{j}} δ_{il} \frac{π̃_{l}Y_{l}}{1 + ∑_{k=1}^{n} \frac{ξ_{ik}δ_{ik}π̃_{k}Y_{k}}{1 + ∑_{j=1}^{m} ξ̃_{jk}δ_{jk}π_{j}X_{j}}} .

Multiplying the factors out, we obtain

C_{il} ≈ \frac{π_{i}X_{i}δ_{il}π̃_{l}Y_{l}}{1 + ∑_{j=1}^{m} ξ̃_{jl}δ_{jl}π_{j}X_{j} + ∑_{k=1}^{n} ξ_{ik}δ_{ik}π̃_{k}Y_{k} \frac{1 + ∑_{j=1}^{m} ξ̃_{jl}δ_{jl}π_{j}X_{j}}{1 + ∑_{j=1}^{m} ξ̃_{jk}δ_{jk}π_{j}X_{j}}} .  (2.11)

Notice that this formula coincides with the symmetric formula (2.10), if there is only one predator species. For one prey species, this formula is the same as formula (13) which De Boer and Perelson (1995) obtain for T cell proliferation by a Padé approximation. While the symmetric situation in T cell proliferation may not allow a generalization of the one-antigen many-T-cell-clones scenario to a scenario with many antigen strains and many T cell clones, the asymmetry in the predator-prey relationship makes such a generalization possible. The same holds for mating if one sex vastly outnumbers the other.

Under proportionate mixing (Section 3), we will find a flip-flop approximation which is an estimate from above for C_{il} and is sometimes closer to C_{il} than the approximation (2.10) and which can formally be generalized as

C_{il} ≈ \frac{δ_{il}π_{i}π̃_{l}X_{i}Y_{l}}{1 + max \{ ∑_{k=1}^{n} ξ_{ik}δ_{ik}π̃_{k}Y_{k}, ∑_{j=1}^{m} ξ̃_{jl}δ_{jl}π_{j}X_{j}\}}.  (2.12)

In general, we see that the right hand side of (2.12) is larger than the right hand side of (2.10), but we could not determine under which conditions it is a better approximation. The following compound parameters appear in the last three formulas and the subsequent sections, so we recall their definition in terms of original parameters (cf. (2.4) and (2.7)),

ξ_{jk}δ_{jk} = \frac{γ_{jk}}{ρ_{jk}} + κ\frac{p_{jk}γ_{jk}}{σ_{jk}} ,  ξ̃_{jk}δ_{jk} = \frac{γ_{jk}}{ρ_{jk}} + \frac{p_{jk}γ_{jk}}{σ_{jk}} .  (2.13)

Again, under the condition that the γ_{il} are very small compared with the ρ_{il} and σ_{il} (e.g., the time it takes to locate a potential mate is large compared with the time spent on courtship and the actual mating), we also have global asymptotic stability of the equilibrium (Appendices C and D).

Theorem 2.3 The unique equilibrium of system (2.1), (2.2) is globally asymptotically stable for non-negative solutions (in the sense of Theorem 2.1) if

θ ∑_{l=1}^{n} π̃_{l}Y_{l} + θ̃∑_{i=1}^{m} π_{i}X_{i} < 1,

θ = max_{ik} (\frac{γ_{ik}}{ρ_{ik}} + κ\frac{p_{ik}γ_{ik}}{σ_{ik}} ),

θ̃ = max_{ik} (\frac{γ_{ik}}{ρ_{ik}} + \frac{p_{ik}γ_{ik}}{σ_{ik}} ).

For local asymptotic stability alone, the assumption in Theorem 2.3 can be slightly relaxed.

3 Proportionate mixing

The term proportionate mixing comes from the sexually transmitted disease literature where it means that, apart from different activity levels, the different groups mix homogeneously. By abuse of language, we speak about proportionate mixing in the context of complex formation if the proportions

\frac{γ_{jl}}{ρ_{jl}} , \frac{p_{jl}γ_{jl}}{σ_{jl}}

are independent of j and l. In other words, if there exist constants c, c̃ > 0 such that

γ_{jl} = cρ_{jl} = c̃\frac{σ_{jl}}{p_{jl}} .

This means, in the context of predator prey interaction, e.g., that there is a species independent proportion between the time it takes the predator to find the prey and the time spent on trying to kill the prey and the time spent on eating and digesting the prey multiplied with the probability of killing the prey without being killed.

By (2.11), proportionate mixing implies that the compound parameters ξ_{jk}δ_{jk} and ξ̃_{jk}δ_{jk} do not depend on j and k. For the subsequent considerations it is actually sufficient that ξ_{jk}δ_{jk} do not depend on j, while ξ̃_{jk}δ_{jk} do not depend on k. It then follows from (2.6) that u_{i} is independent of i, w_{l} is independent of l and K_{il} = K is independent of i and l. Further K satisfies

K = (1 − ȳK)(1 − x̄K)

with

x̄ = ∑_{j=1}^{m} π_{j}X_{j}ξ̃_{jl}δ_{jl}, ȳ = ∑_{k=1}^{n} π̃_{k}Y_{k}ξ_{ik}δ_{ik}.  (3.1)

Solving this quadratic equation we obtain

K = \frac{1}{2x̄ȳ} [1 + x̄ + ȳ − \sqrt{(1 + x̄ + ȳ)^{2} − 4x̄ȳ}] .

The following form is more useful,

K = \frac{2}{1 + x̄ + ȳ + \sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}}} .  (3.2)

Among other things it reveals why we have taken this and not the other root which would be infinity for x̄ = ȳ = 0.

Recalling (2.3) we have

D_{il}^{∗} = \frac{p_{il}ρ_{il}}{σ_{il}} C_{il}^{∗}, (3.3)

C_{il}^{∗} = \frac{2δ_{il}π_{i}X_{i}π̃_{l}Y_{l}}{1 + x̄ + ȳ + \sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}}} ,

with x̄ and ȳ from (3.1). We call this formula for C_{il} and D_{il} the complex formation functional response. Apparently it was first derived by Hsu and Fredrickson (1975) in the context of age-dependent mating (so the heterogeneity of the problem came from age structure rather than different subgroups or species) and then rediscovered for predator prey models by Huisman and De Boer (1997) for the case of one prey and one predator species. Heesterbeek and Metz (1993) found a related formula for sexual disease transmission in homosexual populations. Ruxton, Gurney, and de Roos (1992) derived a different predator prey functional relationship with a square root expression; their derivation is also based on Michaelis-Menten type time scale arguments but does not really involve the formation of prey predator complexes but aggressive encounters between different predators.

While formula (3.3) may look complicated, it is still elegant because of its symmetries and also very informative, because it will give us some indication under which circumstances simpler functional responses are good approximations of the complex formation approach. In the next section, in a predator prey model, we will even see that we may loose interesting qualitative behavior of the model dynamics when we replace the complex formation functional response by an approximation.

4 Estimates and comparisons

The formulas (3.2) and (3.3) in the previous section show that the function

ζ(x̄, ȳ) = 1 + x̄ + ȳ + \sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}}

plays a central role in the complex formation functional response. In this section we investigate whether it can be approximated by a simpler expression. The following estimates hold.

Lemma 4.1 2(1 + max{x̄, ȳ}) ≤ ζ(x̄, ȳ) ≤ 2(1 + x̄ + ȳ).

Proof: The second estimate is obvious. As for the first,

ζ(x̄, ȳ) ≥ 1 + x̄ + ȳ + \sqrt{1 + 2|x̄ − ȳ| + (x̄ − ȳ)^{2}} = 1 + x̄ + ȳ + 1 + |x̄ − ȳ|.

Both estimates are candidates for approximations. Choosing the upper estimate would lead to the functional response suggested by Beddington (1975) and DeAngelis et al. (1975). In order to find out under which conditions one estimate is better than the other, we derive the following estimates for their distances from ζ.

Lemma 4.2

\frac{2x̄ȳ}{1 + x̄ + ȳ} ≤ 2(1 + x̄ + ȳ) − ζ(x̄, ȳ) ≤ \frac{2x̄ȳ}{1 + max\{x̄, ȳ\}}, (a)

\frac{2 min\{x̄, ȳ\}}{1 + max\{x̄, ȳ\}} ≤ ζ(x̄, ȳ) − 2(1 + max\{x̄, ȳ\}) ≤ \frac{2 min\{x̄, ȳ\}}{1 + |x̄ − ȳ|}. (b)

Proof: Notice that

2(1 + x̄ + ȳ) − ζ(x̄, ȳ) = 1 + x̄ + ȳ − \sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}}
= \frac{4x̄ȳ}{1 + x̄ + ȳ + \sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}}} .

Replacing x̄ + ȳ by |x̄ − ȳ| or vice versa yields the two estimates in (a). Further

ζ(x̄, ȳ) − 2(1 + max\{x̄, ȳ\}) = \sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}} − (1 + |x̄ − ȳ|)
= \frac{2(x̄ + ȳ − |x̄ − ȳ|)}{\sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}} + (1 + |x̄ − ȳ|)}
= \frac{4 min\{x̄, ȳ\}}{\sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}} + (1 + |x̄ − ȳ|)}

Replacing x̄ + ȳ by |x − y| or vice versa under the square root yields the two estimates in (b).

Approximations

Lemma 3.2 (a) shows that, if x̄ and ȳ are small, ζ(x̄, ȳ) ≈ 2(1 + x̄ + ȳ) and we can approximate C by a functional response à la Beddington and DeAngelis et al.

C_{il} ≈ \frac{δ_{il}π_{i}X_{i}π̃_{l}Y_{l}}{1 + x̄ + ȳ} .  (3.4)

The functional response (3.4) can also be derived as an approximation of the square root formula obtained by Ruxton et al. (1992) and from effects of spatial grouping (Cosner et al., 1999).

However, if one of x̄ and ȳ is small and the other large, the following approximation of flip-flop type is better than the approximation (3.4),

C_{il} ≈ \frac{δ_{il}π_{i}X_{i}π̃_{l}Y_{l}}{1 + max\{x̄, ȳ\}}.  (3.5)

As seen from Lemma 3.2, the flip-flop response is no good approximation at all, if x̄ and ȳ are about equal, even if both are small. As we will see, the approximation (3.4) has the advantage to extend to general mixing (see (2.8) and Appendix B), while this is not clear for the flip-flop approximation. In predatorprey models, if ȳ < x̄ (i.e., the number of predators is small compared to the number of prey), the flip-flop approximation reduces to the Michaelis-Menten (1913) alias Monod (1942, 1950) alias Holling type 2 (1965, 1966) functional response.

Compared with the functional response (3.4) (and the complex formation functional response) the Michaelis-Menten response neglects direct or indirect interaction of predators with each other. In Ruxton et al. (1992) this interaction is direct and originates from aggressive encounters between predators, while in the complex formation approach the interaction is indirect as a prey engaged by one predator is assumed to be not available for others. See also the comparison in Heesterbeek and Metz (1993).

If we are dealing with one species on each side only, we have

C = \frac{2}{δξξ̃}\frac{x̄ȳ}{ζ(x̄, ȳ)}.

If one of x̄ and ȳ is large and the other small, the flip-flop approximation (3.5) can be further approximated by

C ≈ \frac{1}{δξξ̃}\frac{x̄ȳ}{max\{x̄, ȳ\}} = \frac{1}{δξξ̃} min\{x̄, ȳ\}

which is one of the common functional responses in marriage models. This seems to be the only instance where the approximation of the complex formation functional response leads to a ratio-dependent expression; it would be rather artificial, though, to rewrite the minimum function in ratio-dependent form. There does not seem to be a justification for the harmonic mean (another common marriage function) in this context. In order to obtain it from the functional response (3.4),

C ≈ \frac{1}{δξξ̃}\frac{x̄ȳ}{1 + x̄ + ȳ},

we would need to assume that x̄ + ȳ is large, while in justifying (3.4) we needed to assume that this sum is small.

5 A prey model with a generalist predator

This section presents a case study as to what extent the complex formation functional response may be replaced by the response à la Beddington (1975) and DeAngelis et al. (1975) or by the flip-flop response without changing the qualitative behavior of a model.

In a prey model with a generalist predator it is assumed that the predator has enough alternative prey that it does not profit from eating this specific prey. We also assume that the prey never kills the predator. If the prey species, whose size is denoted by X, obeys a standard logistic equation in the absence of the predator, the model reads

X′ = r̃X − ν̃X^{2} − pρC

where C denotes the first stage of the prey-predator complex. Recall that ρC is the exit rate of the first stage and pρC is the rate at which prey is killed. p is the probability at which individual prey is killed at the end of the first stage. By (3.3) and (3.1),

C = \frac{2δπXπ̃Y}{1 + x̄ + ȳ + \sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}}} .

x̄ = πXξ̃δ , ȳ = π̃Y ξδ .

One may wonder whether r̃X should be replaced by r̃(X − C) and ν̃X^{2} by ν̃(X −C)^{2} because prey in a complex would not give birth to surviving offspring or take part in overcrowding effects (cf. Huisman and De Boer, 1997). However, the quasi-steady-state approach is only justified under the assumption that the predation process is much faster than the demographic processes in the prey population which means that the C terms arising from this modification can also be neglected in comparison with pρC.

Introducing new parameters and a new time leads to the following equation for x̄:

\frac{x̄′}{x̄} = r − νx̄ − \frac{2}{1 + x̄ + ȳ + \sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}}} .  (4.1)

Notice that ȳ which is proportional to the size of the predator species, Y, is just a, though crucial, parameter.

We want to compare the qualitative behavior of the solutions of this equation with the one of solutions to the approximating equations

\frac{x̄′}{x̄} = r − νx̄ − \frac{1}{1 + x̄ + ȳ},  (4.2)

and

\frac{x̄′}{x̄} = r − νx̄ − \frac{1}{1 + max\{x̄, ȳ\}}.  (4.3)

Equilibria of equation (4.1) are intersection points of the decreasing line (1/2)(r− νx̄) and the graph of the function

ϕ(x̄) = \frac{1}{1 + x̄ + ȳ + \sqrt{1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2}}} .

One easily sees that ϕ is strictly decreasing. As far as convexity is concerned, there are two cases (Appendix E):

Case 1: ȳ ≤ 1. Then ϕ is strictly convex.

Case 2: ȳ > 1. There exists x̄^{♯} > 0 such that ϕ is strictly concave on [0, x̄^{♯}] and strictly convex on [x̄^{♯}, ∞).

In the second case, we could not determine x̄^{♯} precisely, but at least locate it in the interval (ȳ− 1, (5/3)ȳ− 1). In the first case, the functional response in (4.2) has the same qualitative features as the complex formation response in (4.1), while in the second case the flip-flop response in (4.3) caricatures, though in an exaggerated way, the features of the complex formation response.

In particular, the square root response can display the following scenario in case 2 which cannot be produced by the response in (4.2) or any of the traditional models for generalist predators; it can be produced, though, by the flip-flop response in (4.3):

Our scalar predator prey model can have as many as three non-trivial equilibria 0 < X_{1}^{∗} < X_{2}^{∗} < X_{3}^{∗}, associated with the following large-time behavior: Prey populations starting with sizes between 0 and X_{2}^{∗} are closely controlled by the predator and converge towards X_{1}. Prey populations starting with sizes larger than X_{2}^{∗} escape the close control of the predator and converge towards X_{3}^{∗}. There is some impact by the predator in so far as X_{3}^{∗} is lower than the carrying capacity of the prey population in absence of the predator.

In all traditional prey predator models, a generalist predator that closely controls the prey necessarily drives it into extinction, unless the generalist predator lowers its search intensity at low prey densities as known from the Holling type 3 functional response.

6 Sexually transmitted diseases in heterosexual populations

Many attempts have been made to model the mixing behavior in multi-group populations in models for sexually transmitted diseases (e.g., Castillo-Chavez, Busenberg, 1991, Castillo-Chavez et al., 1994, Jacquez et al., 1995, and the literature mentioned there). The complex formation approach offers the possibility to model the functional form of contacts between the various groups from first principles. In this section we illustrate how complex formation can be built into epidemic models for heterosexual multi-group populations and how the basic reproductive number, R_{0}, can be determined in the case of separable mixing in spite of the seemingly overwhelming complexity of the model.

Disease transmission occurs, if at all, during the consumption (mating) stage of the complex provided the pair consists of one susceptible and one infective individual. If an infection occurs, the newly infected individual is effectively released into the infective part of its group when the consumption stage ends. So the incidence (rate at which newly infected individuals are introduced into the infective part of their group) in female group i caused by contacts with male group k, is

J_{ik} = \frac{S_{i}}{X_{i}} σ_{ik}q_{ik}D_{ik} \frac{Ĩ_{k}}{Y_{k}} .  (5.1)

Here S_{i} is the number of susceptibles in female group i and \frac{S_{i}}{X_{i}} is the probability that the female group i individual involved in the complex is actually susceptible. Ĩ_{k} is the number of infectives in male group k and \frac{Ĩ_{k}}{Y_{k}} is the probability that the male group k individual involved in the complex is actually infective. q_{ik} is the probability that a contact between a susceptible female from group i and an infective male from group k results in transmission of the disease. By (2.3),

J_{ik} = \frac{S_{i}}{X_{i}} p_{ik}ρ_{ik}q_{ik}C_{ik} \frac{Ĩ_{k}}{Y_{k}} ,

and by (2.5),

J_{ik} = π_{i}S_{i}p_{ik}γ_{ik}q_{ik}K_{ik}π̃_{k}Ĩ_{k},  (5.2)

with the matrix K being the solution of system (2.6).

Similarly, the incidence in male group k caused by contacts with female group i, is given by

J̃_{ik} = \frac{I_{i}}{X_{i}} σ_{ik}q_{ik}D_{ik} \frac{S̃_{k}}{Y_{k}} = π_{i}I_{i}p_{ik}γ_{ik}q̃_{ik}K_{ik}π̃_{k}S̃_{k}.

The incidence in female group i is now given by ∑_{k=1}^{n} J_{ik}, while the incidence in male group l is given by ∑_{j=1}^{m} J̃_{jl}. These incidences can now be used as building blocks in the formulation of epidemic models. Recall that we have an explicit expression for K_{ik} = K in the case of proportionate mixing, (3.1) and (3.2), and the approximate formulas corresponding to (2.8) and (2.10),

K_{il} ≈ \frac{1}{1 + ∑_{k=1}^{n} ξ_{ik}δ_{ik}π̃_{k}Y_{k}} \frac{1}{1 + ∑_{j=1}^{m} ξ̃_{jl}δ_{jl}π_{j}X_{j}} ,

K_{il} ≈ \frac{1}{1 + ∑_{k=1}^{n} ξ_{ik}δ_{ik}π̃_{k}Y_{k} + ∑_{j=1}^{m} ξ̃_{jl}δ_{jl}π_{j}X_{j}} ,

if the sums in the denominator are small. Presumably the second approximation is closer to K_{il} than the first (Appendix B).

The basic reproduction ratio

Let τ_{i} the average length of the infective period for female group i and τ̃_{k} the average length of the infective period for male group k.

Let R_{ik} be the basic reproduction number from male group k to female group i, i.e., the average number of secondary cases one infective male individual in group k can produce in female group i during the whole infective period if the female group i is completely disease-free. (Cf. Jacquez, Simon, Koopman, 1995, where R_{ik} is called the basic reproduction number for the female group i contacts of an infected male in group k.) We have

R_{ik} = J_{ik} \frac{τ̃_{k}}{I_{k}} ,

with X_{i} replacing S_{i} in (5.2), i.e.,

R_{ik} = π_{i}X_{i}p_{ik}γ_{ik}q_{ik}K_{ik}π̃_{k}τ̃_{k}.

Similarly we define the basic reproduction ratio from female group i to male group k, R̃_{ki}, and obtain

R̃_{ki} = π_{i}τ_{i}p_{ik}γ_{ik}q̃_{ik}K_{ik}π̃_{k}Y_{k}.

Let R be the m × n matrix with entries R_{ik} and R̃ the n × m matrix with entries R̃_{ki}. Counting the female groups first, the basic reproduction matrix (or next generation matrix, Diekmann, Heesterbeek, Metz, 1990, 1995) of the disease in the whole population is given by

R^{⋄} = ( \begin{matrix}0 & R \\ R̃ & 0 \end{matrix} ) .

The basic reproduction ratio of the disease, R_{0}, is mathematically defined as the spectral radius of the basic reproduction matrix, R^{⋄}. Epidemiologically, it is the number of secondary cases an infective individual can produce during the whole infectious period if introduced into an otherwise disease-free population (Diekmann, Heesterbeek, Metz, 1990, 1995). R_{0} is associated with a threshold phenomenon for the disease-free state of the population. If R_{0} < 1, the diseasefree state is locally asymptotically stable, i.e., if only a few number of infectives are introduced, the disease dies out. Often, but not always, the disease-free state is even globally asymptotically stable, which means that the disease finally dies out also when many infectives are introduced. If R_{0} > 1, the disease-free state is unstable and often the disease becomes endemic in the population.

We have

(R^{⋄})^{2} = ( \begin{matrix}RR̃ & 0 \\ 0 & R̃R \end{matrix} )

which implies that R_{0}^{2} is the spectral radius of the matrix RR̃. (Cf. Diekmann, Dietz, Heesterbeek, 1991.)

The case of separable mixing

Here we assume that the matrices p_{ik}, γ_{ik}, q_{ik}, q̃_{ik} are separable, i.e., they are of the form α_{i}β_{k}. Since the matrix K_{ik} is separable as well by (2.6), the matrices R and R̃ are also separable, i.e., of the form

R_{jk} = α_{j}β_{k}
R̃_{kj} = α̃_{j}β̃_{k} .

The product of these two separable matrices is also separable and its spectral radius is given by

( ∑_{j=1}^{m} α_{j}α̃_{j})( ∑_{k=1}^{n} β_{k}β̃_{k}).

Hence

R_{0}^{2} = ∑_{j=1}^{m} ∑_{k=1}^{n} R_{jk}R̃_{kj}.

Substituting the expressions above,

R_{0}^{2} = ∑_{j=1}^{m} ∑_{k=1}^{n}(π_{j})^{2}(π̃_{k})^{2}p_{jk}^{2}γ_{jk}^{2}q_{jk}q̃_{jk}K_{jk}^{2}τ_{j}τ̃_{k}X_{j}Y_{k} .

To evaluate K_{jk} we can use appropriate approximations. However, the fixed point theory outlined in Appendix B allows to study K directly as a continuous function of X = (X_{1}, …, X_{m}), Y = Y_{1}, …, Y_{m}) and offers an iterative scheme to calculate convergent numerical approximations for K_{jk} which, at the same time, are estimates from below and above respectively.

7 Discussion

Over the years, the mass action type functional relationship between prey and predators suggested by Lotka and Volterra has been modified in various directions, both by intuitive reasoning (Beddington, 1975, DeAngelis et al., 1975, Arditi and Ginzburg, 1989, and by derivation from first principles. The original derivation by Holling’s (1965, 1966) concentrates on the functional response of one predator to the amount of available prey and does not consider direct or indirect interference of predators. Direct interference of predator can be built into the original Holling model by assuming that the time available for searching prey is shortened by aggressive encounters with other predators ( Ruxton et al., 1992) or that predators must divide the prey among themselves ( Cosner et al., 1999). The complex formation approach discussed here (and by Borghans et al., 1996, De Boer and Perelson, 1995, and Huisman and De Boer for predatorprey models and/or T cell production) assumes an indirect interference between predators in so far as prey that is engaged or handled by a predator is off limits for other predators. This does not necessarily exclude pack hunting or group defense; if the pack or group sizes are rather constant, one can apply the complex formation argument to packs or groups instead of single individuals. The complex formation approach relies on a time scale argument which justifies taking quasi-steady-states. If there is only one prey and one predator species or if the species mix homogeneously, the quasi-steady-state equations can be solved explicitly and yield similar expressions to those obtained by Ruxton et al. (1992). Hsu and Fredrickson (1975) seem to be the first to find this explicit solution. They derive it from first principles in the context of pair formation and mating, while most other pair formation functions are determined phenomenologically to fit certain desired features (see Hadeler et al., 1988, for a discussion and the history of various marriage functions). It has been noticed before (De Boer, Perelson, 1995, Borghans et al., 1996, Huisman, De Boer, 1997) that the functional relationship suggested by Beddington (1975) and DeAngelis et al. (1975), c\frac{x̄ȳ}{1+x̄+ȳ}, with x̄, ȳ being scaled amounts of predator and prey, is often is a reasonable approximation of the quasi-steady-state of complex formation. This functional relationship can also be obtained by considering space-limited predation ( Cosner et al., 1999). However, if there are few prey and many predators or many prey and few predators, the flip-flop response c\frac{x̄ȳ}{1+max\{x̄,ȳ\}} is a better approximation of the quasi-steady state in complex formation (Section 3). If max{x̄, ȳ} is large compared to one, the flip-flop relationship takes the form of a familiar marriage function, c min{x̄, ȳ}. The complex formation approach does not provide ratio-dependent functional relationships, except in this special case which can artificially be rewritten in a ratio-dependent form.

While the above-mentioned approximations are appropriate in certain domains of the functional relationship, they are not so suitable in others. To illustrate this we consider the impact of a generalist predator on a prey population and show that with either approximation we do not obtain the full spectrum of quantitative behavior of the complex formation model (Section 4). This is a message to keep in mind when proceeding to multi-species or multi-group models where an explicit formula for the quasi-steady state is not available.

The continuing spread of old and emergence of new sexually transmitted diseases have made it necessary to model the contact structure in very heterogeneous populations. While it is possible to design the relationship of sexual contacts between various groups as a function of group sizes along certain consistency requirements and conservation laws (e.g., Castillo-Chavez, Busenberg, 1991, and Castillo-Chavez et al., 1994), the complex formation approach opens a very natural and attractive modeling avenue to a derivation from first principles (Section 5). It faces the shortcoming, however, that the quasi-steady-states cannot be determined explicitly for a multiple group situation except for homogeneous or a special form of proportionate mixing. While this difficulty was already noticed in the pair formation and predator prey context, De Boer and Perelson (1995) and Huisman and De Boer (1997) were able to present approximate solutions in case that there is one prey species and many predator species or one predator and many prey species. We extend these formulas to the case of many prey and predator species (Section 2 and Appendix B). Due to the inherent asymmetry in a prey-predator relationship, they are asymmetric generalizations of the functional forms found by Beddington (1975) and DeAngelis et al. (1975).

For sexually transmitted diseases, in heterosexual populations with multiple groups of men and women, we suggest an approximation which is symmetric in men and women. Miraculously in view of the complexity of the model, in the case of separable mixing one can derive a formula for the basic reproduction number of the disease in terms of the implicit solutions of the quasi-steady-state equations which offers a very convenient way to determine it precisely by straightforward numerical iterations (Section 5).

Acknowledgements The first author thanks Maia Martcheva for drawing his attention to the paper by Hsu and Fredrickson (1975) which had apparently been overlooked in the more recent literature. For further bibliographic hints he thanks Chris Cosner and Hans (J.A.J.) Metz.

Appendix A: Global existence and uniqueness of solutions

Substituting (2.1) into (2.2) yields a system of ordinary differential equations in the space of real m × n matrices which can be identified with ℝ^{nm}. Since the vector field is everywhere defined and continuously differentiable, standard ODE theory (e.g. Hale, 1980, Section I.1 to I.3) implies that, for all initial data, there exists a local solution on an interval [0, b) which is maximal in the sense that the norm of the solution becomes unbounded as t → b in case that b < ∞.

Define a closed convex set C as the set of those C_{il} ≥ 0 such that x_{i}, y_{l} ≥ 0 with these being given by (2.1). From (2.1) we obtain that, on [0, b),

x′_{i} = − ∑_{k=1}^{n} C′_{ik} − κ ∑_{k=1}^{n} D′_{ik}

= ∑_{k=1}^{n} ρ_{ik}[1 − κp_{ik}]C_{ik} + κ ∑_{k=1}^{n} σ_{ik}D_{ik} − ∑_{k=1}^{n} γ_{ik}x_{i}y_{k},

y′_{l} = ∑_{j=1}^{m} ρ_{jl}[1 − p_{jl}]C_{jl} + ∑_{j=1}^{m} σ_{jl}D_{jl} − ∑_{j=1}^{m} γ_{jl}x_{j}y_{l} .

It follows from Proposition B.7 in Smith, Waltman (1995) that the convex set C is forward invariant under the solution. By (2.1), solutions C_{ik} ≥ 0 with x_{i}, y_{l} ≥ 0 are bounded, which implies that solutions which satisfy these constraints initially are defined for all forward times and satisfy the constraints forever.

Appendix B: Uniqueness of the steady state

We will show that the system (2.6) for given X = (X_{1}, …, X_{m}), Y = (Y_{1}, …, Y_{n}), has a unique solution such that K_{il} and the two factors on the right hand side are non-negative. From (2.6),

K_{il} = u_{i}w_{l}

u_{i} = 1 − ∑_{k=1}^{n} α_{ik}K_{ik} (B1)

w_{l} = 1 − ∑_{j=1}^{m} α̃_{jl}K_{jl},

where we are interested in u_{i} and w_{l} being non-negative. Here

α_{il} = ξ_{il}δ_{il}π̃_{l}Y_{l} , α̃_{il} = ξ̃_{il}δ_{il}π_{i}X_{i} .  (B2)

In (B1) we substitute the first into the second and third equation,

u_{i} = 1 − u_{i}( ∑_{k=1}^{n} α_{ik}w_{k}) , w_{l} = 1 − w_{l}( ∑_{j=1}^{m} α̃_{jl}u_{j}) .

We solve for u_{i} and w_{l},

u_{i} = \frac{1}{1 + ∑_{k=1}^{n} α_{ik}w_{k}} , w_{l} = \frac{1}{1 + ∑_{j=1}^{m} α̃_{jl}u_{j}} .  (B3)

Notice that any non-negative solution u = (u_{1}, …, u_{m}), w = (w_{1}, …, w_{n}) of (B3) satisfies

u_{i} ≤ 1, w_{l} ≤ 1,

and that K_{il} = u_{i}w_{l} provides a solution to (2.6).

To investigate whether the system (B3) has solutions, we substitute the second part of the system into the first:

u_{i} = \frac{1}{1 + ∑_{k=1}^{n} \frac{α_{ik}}{1 + ∑_{j=1}^{m} α̃_{jk}u_{j}}} =: G_{i}(u).  (B4)

Let G : [0, ∞)^{m} → [0, ∞)^{m} be the mapping formed by taking the G_{i} as components. Apparently it is sufficient to find a fixed point of G.

We introduce an order on ℝ^{m} by setting

x ≤ y ⟺ x_{i} ≤ y_{i}, i = 1, …, m,

for vectors x = (x_{1}, …, x_{m}), y = (y_{1}, …, y_{m}). We see that G is an increasing map, i.e.,

G(u) ≤ G(v) whenever u ≤ v.

We also notice that G(u) ≤ 1 for all u ∈ [0, ∞)^{m} where 1 = (1, …, 1). Let [0, 1] denote the generalized order interval of those vectors u with 0 ≤ u ≤ 1. G maps the closed convex bounded set [0, 1] into itself and is continuous. Hence G has a fixed point u in [0, 1] by Brouwer’s fixed point theorem.

Actually G maps [0, ∞)^{m} into the order interval [G(0), 1], and G(0) ≥ δ1 for some δ > 0. In order to show that there exists exactly one fixed point of G in [0, ∞)^{m}, we use the concavity (or sub-homogeneity) technique of Krasnosel’skii (1964), Section 6.1.

Let v be a positive vector and ξ ∈ (0, 1). Then

G_{i}(ξv) = \frac{1}{1 + ∑_{k=1}^{n} α_{ik}(\frac{1}{1+∑_{j=1}^{m} α̃_{jk}ξv_{j}} )} ≥ \frac{1}{1 + ∑_{k=1}^{n} α_{ik}(\frac{1}{ξ+∑_{j=1}^{m} α̃_{jk}ξv_{j}} )}

= \frac{1}{1 + \frac{1}{ξ} ∑_{k=1}^{n} α_{ik}(\frac{1}{1+∑_{j=1}^{m} α̃_{jk}v_{j}} )} = \frac{ξ}{ξ + ∑_{k=1}^{n} α_{ik}(\frac{1}{1+∑_{j=1}^{m} α̃_{jk}v_{j}} )}

> ξ \frac{1}{1 + ∑_{k=1}^{n} α_{ik}(\frac{1}{1+∑_{j=1}^{m} α̃_{jk}v_{j}} )} = ξG_{i}(v).

This implies that G is 1-concave in the terminology of Krasnosel’skii and so the uniqueness of fixed points and the convergence of successive approximations (Krasnosel’skii, 1964, 6.1.3 and 6.1.7).

In particular the following monotone sequences converge towards the unique fixed point,

u(0) = 0, u(j + 1) = G(u(j)),
v(0) = 1, v(j + 1) = G(v(j)) .

The first sequence is monotone increasing and the second monotone decreasing, they provide improving upper and lower estimates for the fixed point. Estimate (2.8) and the subsequent estimate are obtained from G(0) and G(1) and analogous considerations for w in (B3) rather than u (see also the next paragraph).

A Ross type estimate and approximation

From (B3) we obtain an estimate from below for K_{il},

K_{il} = u_{i}w_{l} ≥ \frac{1}{1 + ∑_{k=1}^{n} α_{ik}} \frac{1}{1 + ∑_{j=1}^{m} α̃_{jl}}  (B5)

which leads to estimate (2.8) involving an expression reminiscent of Ross solutions (Castillo-Chavez, Busenberg, 1991). If the α’s and α̃’s are small, (B5) actually provides a reasonable approximation of K_{il}.

A Beddington type approximation

Dropping the quadratic terms in the denominator in (B5) we obtain the following approximation which is an extension of the one suggested by Beddington (1975) and DeAngelis et al. (1975),

K_{il} ≈ \frac{1}{1 + ∑_{k=1}^{n} α_{ik} + ∑_{j=1}^{m} α̃_{jl}} .  (B6)

The following consideration shows that the approximation (B6) may be even better than the Ross type approximation. By repeated use of (B1),

\frac{1}{K_{il}} = \frac{1}{u_{i}w_{l}} = (1 + ∑_{k=1}^{n} α_{ik}w_{k}) (1 + ∑_{j=1}^{m} α̃_{jl}u_{j})

= 1 + ∑_{k=1}^{n} α_{ik}w_{k} + ∑_{j=1}^{m} α̃_{jl}u_{j} + ∑_{k=1}^{n} ∑_{j=1}^{m} α_{ik}w_{k}α̃_{jl}u_{j}

= 1 + ∑_{k=1}^{n} α_{ik}(1 − w_{k} ∑_{j=1}^{m} α̃_{jk}u_{j}) + ∑_{j=1}^{m} α̃_{jl}(1 − u_{j} ∑_{k=1}^{n} α_{jk}w_{k})

+ ∑_{k=1}^{n} ∑_{j=1}^{m} α_{ik}w_{k}α̃_{jl}u_{j}.

Hence

K_{il} = \frac{1}{1 + ∑_{k=1}^{n} α_{ik} + ∑_{j=1}^{m} α̃_{jl} + ∑_{k=1}^{n} ∑_{j=1}^{m} [α_{ik}α̃_{jl} − α_{ik}α̃_{jk} − α_{jk}α̃_{jl}] u_{j}w_{k}} .

If there is only one prey or one predator species, the double sum in the denominator is negative, so the approximation (B6) is also a lower estimate of the complex formation expression, but a better one than the Ross approximation. In the case of proportionate mixing we clearly see that, under all circumstances, the approximation (B6) is better than the Ross type approximation because, by Lemma 3.1,

K ≥ \frac{1}{1 + x̄ + ȳ} ≥ \frac{1}{1 + x̄} \frac{1}{1 + ȳ}.

We conclude by mentioning that, if the γ_{il} are very small compared with the ρ_{il} and σ_{il} (i.e., the time needed to find one prey is large compared with the time needed to eventually hunt it down, kill it, eat and digest it), then the α_{il} and α̃_{il} are small and the products α_{ik}α̃_{jl} are very small.

A specific predator-prey approximation

The asymmetry in a predator-prey relationship (κ = 0 in (2.7), only the predator takes active part in the second stage of the complex) has the consequence that ξ_{ik} may be considerably smaller than ξ̃_{ik}. Further the readiness of the predator to hunt should on the average be lower than the availability of the prey to be hunted so that it is reasonable to assume that α_{il} is considerably smaller than α̃_{il} in (B2). So in (B3) we guess that u_{i} is rather close to 1, while this may not be the case for w_{l}. So, by (B3) and (B4), we choose the approximations

u_{i} ≈ \frac{1}{1 + ∑_{k=1}^{n} α_{ik}(\frac{1}{1+∑_{j=1}^{m} α̃_{jk}} )}

w_{l} ≈ \frac{1}{1 + ∑_{j=1}^{m} α̃_{jl}} .

From (B1), multiplying the factors out,

K_{il} ≈ \frac{1}{1 + ∑_{j=1}^{m} α̃_{jl} + ∑_{k=1}^{n} α_{ik} \frac{1 + ∑_{j=1}^{m} α̃_{jl}}{1 + ∑_{j=1}^{m} α̃_{jk}}} .  (B7)

Notice that this yields the approximation (B6) if there is only one predator species. Going back to the original variables yields (2.11).

Appendix C: Local stability of the steady state

Replacing x_{i}, y_{i}, C_{il}, D_{il} by x_{i} + x_{i}^{∗}, y_{i} + y_{i}^{∗}, C_{il} + C_{il}^{∗}, D_{il} + D_{il}^{∗} and dropping the higher order terms we obtain the linearization of the system (2.1), (2.2) at the steady state,

0 = x_{i} + ∑_{k=1}^{n}[C_{ik} + κD_{ik}]

0 = y_{l} + ∑_{j=1}^{m}[C_{jl} + D_{jl}]

C′_{il} = γ_{il}x_{i}^{∗}y_{l} + γ_{il}x_{i}y_{l}^{∗} − ρ_{il}C_{il}

D′_{il} = p_{il}ρ_{il}C_{il} − σ_{il}D_{il}.

Substituting the first two equations into the last two, we obtain the following linear ODE system,

C′_{il} = −γ_{il}x_{i}^{∗} ∑_{j=1}^{m}[C_{jl} + D_{jl}] − γ_{il}y_{l}^{∗} ∑_{k=1}^{n}[C_{ik} + κD_{ik}] − ρ_{il}C_{il}

D′_{il} = p_{il}ρ_{il}C_{il} − σ_{il}D_{il}.

The algebraic system for the associated eigenvalues λ and their eigenvectors has the form

λC_{il} = −γ_{il}x_{i}^{∗} ∑_{j=1}^{m}[C_{jl} + D_{jl}] − γ_{il}y_{l}^{∗} ∑_{k=1}^{n}[C_{ik} + κD_{ik}] − ρ_{il}C_{il}

λD_{il} = p_{il}ρ_{il}C_{il} − σ_{il}D_{il}

where at least one C_{il} or D_{il} is different from 0. Usually one proceeds from this finite linear system to the associated characteristic polynomial and studies the location of its zeros, often using the Routh-Hurwitz criterion. Since our system is arbitrarily large and has an unusual form, we follow a different approach used before in Thieme (1985) and Hethcote, Thieme (1985). Solving for D_{il} and reorganizing the terms, we have

\frac{λ + ρ_{il}}{γ_{il}} C_{il} = −x_{i}^{∗} ∑_{j=1}^{m} C_{jl} [1 + \frac{p_{jl}ρ_{jl}}{λ + σ_{jl}} ] − y_{l}^{∗} ∑_{k=1}^{n} C_{ik} [1 + κ\frac{p_{ik}ρ_{ik}}{λ + σ_{ik}} ]

where at least one C_{il} is different from 0. Taking absolute values,

\frac{|λ + ρ_{il}|}{γ_{il}}|C_{il}| ≤ x_{i}^{∗} ∑_{j=1}^{m} |C_{jl}| [1 + \frac{p_{jl}ρ_{jl}}{|λ + σ_{jl}|}] + y_{l}^{∗} ∑_{k=1}^{n} |C_{ik}| [1 + κ\frac{p_{ik}ρ_{ik}}{|λ + σ_{ik}|}] .

The principle of linearized stability implies that the steady state is locally asymptotically stable if all eigenvalues have strictly negative real part. To show the latter we suppose that ℜλ ≥ 0. Then

\frac{ρ_{il}}{γ_{il}} |C_{il}| ≤ x_{i}^{∗} ∑_{j=1}^{m} |C_{jl}| [1 + \frac{p_{jl}ρ_{jl}}{σ_{jl}} ] + y_{l}^{∗} ∑_{k=1}^{n} |C_{ik}| [1 + κ\frac{p_{ik}ρ_{ik}}{σ_{ik}} ] .

Set ξ = max_{il} \frac{|C_{il}|}{C_{il}^{∗}}. Then ξ > 0. Choose i, l such that |C_{il}| = ξC_{il}^{∗}. For such a pair i, l we have

\frac{1}{δ_{il}} ξC_{il}^{∗} ≤ x_{i}^{∗} ∑_{j=1}^{m} ξC_{jl}^{∗} [1 + \frac{p_{jl}ρ_{jl}}{σ_{jl}} ] + y_{l}^{∗} ∑_{k=1}^{n} ξC_{ik}^{∗} [1 + κ\frac{p_{ik}ρ_{ik}}{σ_{ik}} ] .

Dividing by ξ and using that C_{il}^{∗} = δ_{il}x_{i}^{∗}y_{l}^{∗},

x_{i}^{∗}y_{l}^{∗} ≤ x_{i}^{∗} ∑_{j=1}^{m} δ_{jl}x_{j}^{∗}y_{l}^{∗} [1 + \frac{p_{jl}ρ_{jl}}{σ_{jl}} ] + y_{l}^{∗} ∑_{k=1}^{n} δ_{ik}x_{i}^{∗}y_{k}^{∗} [1 + κ\frac{p_{ik}ρ_{ik}}{σ_{ik}} ] .

We divide by x_{i}^{∗} and y_{l}^{∗} and use (2.7),

1 ≤ ∑_{j=1}^{m} δ_{jl}x_{j}^{∗}ξ̃_{jl} + ∑_{k=1}^{n} δ_{ik}y_{k}^{∗}ξ_{ik}.

Using transformation (2.5),

1 ≤ ∑_{j=1}^{m} δ_{jl}π_{j}X_{j}u_{j}ξ̃_{jl} + ∑_{k=1}^{n} δ_{ik}π̃_{k}Y_{k}w_{k}ξ_{ik}.

Since u_{j}, w_{k} < 1, we obtain a contradiction by assuming

∑_{j=1}^{m} δ_{jl}π_{j}X_{j}ξ̃_{jl} + ∑_{k=1}^{n} δ_{ik}π̃_{k}Y_{k}ξ_{ik} ≤ 1 ∀i = 1, …, m, l = 1, …, n.  (C1)

So this is a sufficient condition for local asymptotic stability of the steady state. Retracing the definition of α_{il} and α̃_{il} we find that it is satisfied if the time needed to form a complex is sufficiently long compared with the durations of the two stages of the complex.

Sharper, but more complicated, conditions can be derived using

u_{i} ≤ \frac{1}{1 + ∑_{k=1}^{n} α_{ik}(\frac{1}{1+∑_{j=1}^{m} α̃_{jk}} )}

and a similar estimate for w_{l}.

Appendix D: Global asymptotic stability of the steady state

If f : [0, ∞) → ℝ is bounded, let

f^{∞} = lim sup_{t→∞} f(t), f_{∞} = lim inf_{t→∞} f(t).

By the Fluctuation Lemma (Hirsch et al., 1995) there exists a sequence t_{j} → ∞ (which depends on i and l) such that

C_{il}(t_{j}) → C_{il}^{∞}, C′_{il}(t_{j}) → 0, j → ∞.

The first equation in (2.2) implies the estimate

C_{il}^{∞} ≤ δ_{il}x_{i}^{∞}y_{l}^{∞},  (D1)

with δ_{ij} in (2.4). Applying the same procedure to the second equation in (2.2),

D_{il}^{∞} ≤ \frac{p_{il}ρ_{il}}{σ_{il}} C_{il}^{∞}.  (D2)

Using analogous arguments for the limit inferior we have

C_{il∞} ≥ δ_{il}x_{i∞}y_{l∞} , D_{il∞} ≥ \frac{p_{il}ρ_{il}}{σ_{il}} C_{il∞}.  (D3)

Subtracting the first inequality in (D3) from (D1),

0 ≤ C_{il}^{∞} − C_{il∞} ≤ δ_{il}x_{i}^{∞}(y_{l}^{∞} − y_{l∞}) + δ_{il}(x_{i}^{∞} − x_{i∞})y_{l∞}.

From (2.1),

0 ≤ C_{il}^{∞} − C_{il∞} ≤ δ_{il}π_{i}X_{i}(y_{l}^{∞} − y_{l∞}) + δ_{il}π̃_{l}Y_{l}(x_{i}^{∞} − x_{i∞}).  (D4)

Further, from (2.1),

x_{i}^{∞} ≤ π_{i}X_{i} − ∑_{k=1}^{n}[C_{ik∞} + κD_{ik∞}] ≤ π_{i}X_{i} − ∑_{k=1}^{n} ξ_{ik}C_{ik∞},

y_{l}^{∞} ≤ π̃_{l}Y_{l} − ∑_{j=1}^{m} ξ̃_{jl}C_{jl∞}, x_{i∞} ≥ π_{i}X_{i} − ∑_{k=1}^{n} ξ_{ik}C_{ik}^{∞}, y_{l∞} ≥ π̃_{l}Y_{l} − ∑_{j=1}^{m} ξ̃_{jl}C_{jl}^{∞}.

Substituting these inequalities into (D4),

0 ≤ C_{il}^{∞} − C_{il∞} ≤ δ_{il}π_{i}X_{i} ∑_{j=1}^{m} ξ̃_{jl}[C_{jl}^{∞} − C_{jl∞}] + δ_{il}π̃_{l}Y_{l} ∑_{k=1}^{n} ξ_{ik}[C_{ik}^{∞} − C_{ik∞}].

Set C_{il}^{∞} − C_{il∞} = δ_{il}Δ_{il}, then

0 ≤ Δ_{il} ≤ π̃_{l}Y_{l} ∑_{k=1}^{n} δ_{ik}ξ_{ik}Δ_{ik} + π_{i}X_{i} ∑_{j=1}^{m} δ_{jl}ξ̃_{jl}Δ_{jl}.

Set Δ̄ = ∑_{i=1}^{m} ∑_{l=1}^{n} Δ_{il}, then

0 ≤ Δ̄ ≤ (θ ∑_{l=1}^{n} π̃_{l}Y_{l} + θ̃∑_{i=1}^{m} π_{i}X_{i}) Δ̄,  (D5)

with

θ = max {δ_{ik}ξ_{ik}; i = 1, …, m, k = 1, …, n} ,
θ̃ = max {δ_{ik}ξ̃_{ik}; i = 1, …, m, k = 1, …, n} .

By (D5), if

θ ∑_{l=1}^{n} π̃_{l}Y_{l} + θ̃∑_{i=1}^{m} π_{i}X_{i} < 1,  (D6)

Δ̄ = 0. Since all the summands in Δ̄ are non-negative, Δ_{ik} = 0 for all i, k and so C_{ik}^{∞} − C_{ik∞} = 0. This implies that C_{ik}(t) converges as t → ∞, with the limit being an equilibrium of (2.1, (2.2) which, as we know from Appendix B, is unique. Notice that (D6) implies (C1). (2.7) and (D6) give the assumption in Theorem 2.3.

Appendix E: Proof of cases 1 and 2 in Section 4

Recall

ϕ(x̄) = \frac{1}{1 + x̄ + ȳ + \sqrt{ψ(x̄)}},

(E1)

ψ(x̄) = 1 + 2(x̄ + ȳ) + (x̄ − ȳ)^{2} = (1 + x̄ − ȳ)^{2} + 4ȳ.

Taking the first and second derivative,

ϕ′(x̄) = −\frac{1 + (ψ(x̄))^{−1/2}(1 + x̄ − ȳ)}{(1 + x̄ + ȳ + (ψ(x̄))^{1/2})^{2}},

−ϕ′′(x̄) ≃ [−(ψ(x̄))^{−3/2}(1 + x̄ − ȳ)^{2} + (ψ(x̄))^{−1/2}] (1 + x̄ + ȳ + (ψ(x̄))^{1/2})
−2[1 + (ψ(x̄))^{−1/2}(1 + x̄ − ȳ)]^{2}.

Here ≃ means that the two expressions it connects have the same sign. Multiplying by (ψ(x̄))^{3/2},

−ϕ′′(x̄) ≃ [−(1 + x̄ − ȳ)^{2} + ψ(x̄)] (1 + x̄ + ȳ + (ψ(x̄))^{1/2})
−2(ψ(x̄))^{1/2} [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ]^{2} .

Since −(1 + x̄ − ȳ)^{2} + ψ(x̄) = 4ȳ by (E1), we have

−ϕ′′(x̄) ≃ 2ȳ(1 + x̄ + ȳ + (ψ(x̄))^{1/2})
−(ψ(x̄))^{1/2} [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ]^{2} =: χ(x̄).

It is easy to see that χ(x̄) → −∞ as x̄ → ∞, i.e., ϕ is convex for large x̄. Further

χ(0) = 4ȳ(1 + ȳ) − 4(1 + ȳ) = 4(1 + ȳ)(ȳ − 1),

(E2)

χ(ȳ − 1) = 2ȳ(2ȳ + \sqrt{4ȳ}) − \sqrt{4ȳ}4ȳ = 4ȳ(ȳ − \sqrt{ȳ}),

so χ(0) and χ(ȳ − 1) have the same sign as ȳ − 1. So, if ȳ > 1, ϕ is concave in neighborhoods of 0 and ȳ − 1.

Lemma. ϕ′′ has no zero in [2ȳ−1, ∞) and ϕ is strictly convex on [2ȳ−1, ∞).

Proof: ϕ′′ has a zero x̄ if and only if χ has the same zero x̄, i.e., if and only if

2ȳ(1 + x̄ + ȳ + (ψ(x̄))^{1/2})

= (ψ(x̄))^{1/2} [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ]^{2}
= (ψ(x̄))^{1/2} [(ψ(x̄)) + 2(ψ(x̄))^{1/2}(1 + x̄ − ȳ) + (1 + x̄ − ȳ)^{2}]
= (ψ(x̄))^{1/2}(2(1 + x̄ − ȳ)^{2} + 4ȳ) + 2((1 + x̄ − ȳ)^{2} + 4ȳ)(1 + x̄ − ȳ).

Here we have used (E1) several times. Reorganizing terms we find that ϕ′′(x̄) = 0 if and only if

(ψ(x̄))^{1/2}(2(1 + x̄ − ȳ)^{2} + 2ȳ) = 2ȳ(1 + x̄ + ȳ) − 2((1 + x̄ − ȳ)^{2} + 4ȳ)(1 + x̄ − ȳ).

Simplifying

(ψ(x̄))^{1/2}((1 + x̄ − ȳ)^{2} + ȳ)
= −(1 + x̄ − ȳ)^{3} − ȳ(4(1 + x̄ − ȳ) − (1 + x̄ + ȳ))
= −(1 + x̄ − ȳ)^{3} − ȳ(3(x̄ + 1) − 5ȳ).

The left hand side of this equation is positive. Since the right hand side is negative for x̄ ≥ (5/3)ȳ − 1, there is no equality for x̄ ≥ (5/3)ȳ − 1 and hence no zero of ϕ′′ with x̄ ≥ (5/3)ȳ − 1. Since ϕ′′ is strictly positive for large x̄, we conclude that ϕ′′ is strictly positive on [(5/3)ȳ − 1, ∞).

Since ϕ′′(ȳ− 1) < 0 by (E.2), ϕ′′ must change sign somewhere between ȳ− 1 and (5/3)ȳ− 1. It follows from the next lemma that this is the only sign change of ϕ′′. This finishes the proof of case 2 in Section 4.

Lemma. χ is strictly decreasing on [0, ∞).

Proof: We first show that χ is strictly decreasing on [ȳ − 1, ∞). To this end we consider

χ′(x̄) = 2ȳ(1 + (ψ(x̄))^{−1/2}(1 + x̄ − ȳ))
−(ψ(x̄))^{−1/2}(1 + x̄ − ȳ) [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ]^{2}

−(ψ(x̄))^{1/2}2 [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ] ((ψ(x̄))^{−1/2}(1 + x̄ − ȳ) + 1).

Multiplying by (ψ(x̄))^{1/2} and reorganizing the terms,

χ′(x̄) ≃ 2ȳ((ψ(x̄))^{1/2} + 1 + x̄ − ȳ) − [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ]^{3}
−(ψ(x̄))^{1/2} [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ)]^{2} .

Since

(ψ(x̄))^{1/2} ≥ 1 + |x̄ − ȳ|,
ψ(x̄))^{1/2} + 1 + x̄ − ȳ ≥ 2 + |x̄ − ȳ| + x̄ − ȳ = 2(1 + [x̄ − ȳ]_{+}),

we have

χ′(x̄) ≃ 2ȳ − [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ]^{2} − (ψ(x̄))^{1/2} [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ)] .

(E.3)

We see that χ′(x̄) < 0 if x̄ ≥ ȳ − 1. If ȳ ≤ 1, χ(0) ≤ 0 by (E2) and so that χ(x̄) < 0 for all x̄ > 0. This is case 1 in Section 4.

By (E3), χ′ has a zero x̄ if and only if

2ȳ = [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ]^{2} + (ψ(x̄))^{1/2} [(ψ(x̄))^{1/2} + 1 + x̄ − ȳ)]
= 2ψ(x̄) + 3(ψ(x̄))^{1/2}(1 + x̄ − ȳ) + (1 + x̄ − ȳ)^{2}
= 3(ψ(x̄))^{1/2}(1 + x̄ − ȳ) + 2((1 + x̄ − ȳ)^{2} + 4ȳ) + (1 + x̄ − ȳ)^{2}
= 3(ψ(x̄))^{1/2}(1 + x̄ − ȳ) + 3(1 + x̄ − ȳ)^{2} + 8ȳ.

Reorganizing terms, we see that χ′ has a zero x̄ if and only if

(ψ(x̄))^{1/2}(1 + x̄ − ȳ) = −(1 + x̄ − ȳ)^{2} − 2ȳ.

Squaring both sides we obtain

ψ(x̄)(1 + x̄ − ȳ)^{2} = (1 + x̄ − ȳ)^{4} + 4ȳ(1 + x̄ − ȳ)^{2} + 4ȳ^{2}.  (E.4)

Remembering (E.1), ψ(x̄) = (1 + x̄ − ȳ)^{2} + 4ȳ, the equality (E4) can only hold if ȳ = 0. So χ′ has no zero for ȳ > 0. Since χ′ is strictly negative for large x̄, it is strictly negative everywhere, if ȳ > 0.

References

[1] Arditi, R.; Ginzburg, L.R. (1989): Coupling in predator-prey dynamics: Ratio dependence. J. Theor. Biol. 139, 311-326

[2] Beddington, J.R. (1975): Mutual interference between parasites or predators and its effects on searching efficiency. J. Anim. Ecol. 51, 597-624

[3] Borghans, J.A.M.; De Boer, R.J.; Segel, L.A. (1996): Extending the quasisteady state approximation by changing variables. Bull. Math. Biol. 58, 43-63

[4] Castillo-Chavez, C.; Busenberg, S. (1991): On the solution of the two-sex mixing problem. Differential Equations Models in Biology, Epidemiology and Ecology (S. Busenberg, M. Martelli; eds.), 80-98. Lecture Notes in Biomathematics 92. Springer

[5] Castillo-Chavez, C.; Velasco-Hernandez, J.X.; Fridman, S. (1994): Modeling contact structures in biology. Frontiers in Mathematical Biology (S.A. Levin, ed.), 454-491. Lecture Notes in Biomathematics 100, Springer

[6] Cosner, C.; DeAngelis, D.L.; Ault, J.S.; Olson, D.B. (1999): Effects of spatial grouping on the functional response of predators. Theor. Pop. Biol. 56, 65-75

[7] DeAngelis, D.L.; Goldstein, R.A.; O’Neill, R. V. (1975): A model for trophic interaction. Ecology 56, 881-892 [8] De Boer, R.J.; Perelson, A.S. (1995): Towards a general function describing

T cell proliferation. J. theor. Biol. 175, 567-576

[9] Diekmann, O.; Dietz, K.; Heesterbeek, J.A.P. (1991): The basic reproduction ratio for sexually transmitted diseases: I. Theoretical considerations. Math. Biosci. 107, 325-339

[10] Diekmann, O.; Heesterbeek, J.A.P.; Metz, J.A.J. (1990): On the definition and the computation of the basic reproduction ratio R_{0} in models for infectious diseases in heterogeneous populations. J. Math. Biol. 28, 365-382

[11] Diekmann, O.; Heesterbeek, J.A.P.; Metz, J.A.J. (1995): The legacy of Kermack and McKendrick. Epidemic Models: Their Structure and Relation to Data (D. Mollison, ed.), 95-115. Cambridge Univ. Press

[12] Edelstein-Keshet, L. (1987): Mathematical Models in Biology. McGraw-Hill

[13] Hadeler, K.P.; Waldstätter, R.; Wörz-Busekros, A. (1988): Models for pair formation in bisexual populations. J. Math. Biol. 26, 635-649

[14] Hale, J.K. (1980): Ordinary Differential Equations. Krieger

[15] Heesterbeek, J.A.P.; Metz, J.A.J. (1993): The saturating contact rate in marriage and epidemic models. J. Math. Biol. 31, 529-539

[16] Hethcote, H.W.; Thieme, H.R. (1985): Stability of the endemic equilibrium in epidemic models with subpopulations. Math. Biosciences 75, 205-227

[17] Hirsch, W.M.; Hanisch, H.; Gabriel, J.P. (1985): Differential equations models for some parasitic infections; methods for the study of asymptotic behavior. Comm. Pure Appl. Math. 38, 733-753

[18] Holling, C.S. (1965): The functional response of predators to prey density and its role in mimicry and population regulation. Mem. Ent. Soc. Canada 45, 1-60

[19] Holling, C.S. (1966): The functional response of invertebrate predators to prey density. Mem. Ent. Soc. Canada 48

[20] Hoppensteadt, F. (1974): Asymptotic stability in singular perturbation problems. II: Problems having matched asymptotic expansion solutions. J. Diff. Eqn. 15, 510-521

[21] Hsu, P.-H.; Fredrickson, A.G. (1975): Population-changing processes and the dynamics of sexual populations. Math. Biosci. 26, 55-78

[22] Huisman, G.; De Boer, R.J. (1997): A formal derivation of the ‘Beddington’ functional response. J. theor. Biol. 185, 389-400 [23] Jacquez, J.A.; Simon, C.P.; Koopman, J. (1995): Core groups and the R_{0}s

for subgroups in heterogeneous SIS and SI models. Epidemic Models: Their Structure and Relation to Data (D. Mollison, ed.), 279-301. Cambridge Univ. Press

[24] Krasnosel’skii, M.A. (1964): Positive Solutions of Operator Equations. Noordhoff

[25] Kuang, Y.; Beretta, E. (1998): Global qualitative analysis of a ratiodependent predator-prey system. J. Math. Biol. 36, 389-406

[26] Lin, C.C.; Segel, L.A. (1974, 1988): Mathematics Applied to Deterministic Problems in the Natural Sciences. Siam

[27] Michaelis, L.; Menten, M.I. (1913): Die Kinetik der Invertinwirkung. Biochem. Z. 49, 333-369

[28] Monod, J. (1942): Recherches sur la croissance des cultures bacteriennes. Herman

[29] Monod, J. (1950): La technique de culture continue, théorie et applications. Annales de l’Institute Pasteur 79, 390-401

[30] Ruxton, G.D.; Gurney, W.S.C; de Roos, A.M. (1992): Interference and generation cycles. Theor. Pop. Biol. 42, 235-253

[31] Smith, H.L.; Waltman, P. (1995): The Theory of Chemostat. Dynamics of Microbial Competition. Cambridge Univ. Press

[32] Thieme, H.R. (1985): Local stability in epidemic models for heterogeneous populations. Mathematics in Biology and Medicine (V. Capasso, E. Grosso, S.L. Paveri-Fontana; eds.), 185-189. Lecture Notes in Biomathematics 57. Springer

Horst R. Thieme & Jinling Yang Department of Mathematics, Arizona State University Tempe, AZ 85287-1804, USA e-mail: h.thieme@asu.edu