\documentclass[twoside]{article}
\usepackage{amsfonts} % used for R in Real numbers
\pagestyle{myheadings}
\markboth{ Complex formation approach in modeling predator prey relations}
{ Horst R. Thieme \& Jinling Yang }
\begin{document}
\setcounter{page}{255}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
Nonlinear Differential Equations, \newline
Electron. J. Diff. Eqns., Conf. 05, 2000, pp. 255--283\newline
http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu or ejde.math.unt.edu (login: ftp)}
\vspace{\bigskipamount} \\
%
On the complex formation approach in modeling predator
prey relations, mating, and sexual disease transmission
%
\thanks{ {\em Mathematics Subject Classifications:} 92D25, 92D30, 92D40.
\hfil\break\indent
{\em Key words:} Complex (pair) formation, functional response,
approximations of Beddington, \hfil\break\indent
Michaelis-Menten alias Monod alias Holling 2 and of Ross type,
proportionate mixing, \hfil\break\indent
generalist predator, disease incidence, basic replacement ratio.
\hfil\break\indent
\copyright 2000 Southwest Texas State University.
\hfil\break\indent Published October 25, 2000. \hfil\break\indent
H.R.T. was partially supported by NSF grants DMS-9403884 and DMS-9706787.
} }
\date{}
\author{Horst R. Thieme \& Jinling Yang \\[12pt]
{\em Dedicated to Alan Lazer} \\ {\em on his 60th birthday }}
\maketitle
\begin{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.
\end{abstract}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\section{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).
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: \begin{itemize}
\item present a unified modeling approach to
mating functions, predator-prey relationships and contact rates
in sexually transmitted diseases (Section 2 and Appendix A)
\item 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)
\item use the explicit solution of the steady state
equations in the case of proportionate mixing to compare the
various approximations (Section 3)
\item 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)
\item 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).
\end{itemize}
\section{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 $\pi_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 $\tilde \pi_l$ is the
average probability of an individual from $Y_l$
to be active. These numbers between $0$ and 1 are called the
{\it 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,
$$\displaylines{
\pi_i X_i = x_i + \sum_{k=1}^n [C_{ik} + \kappa D_{ik}],
\quad i =1, \ldots, m,
\cr
\hfill \llap{(2.1)} \cr
\tilde \pi_l Y_l = y_l + \sum_{j=1}^m [C_{jl} + D_{jl} ],
\quad l =1, \ldots, 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
$\kappa = 1$, while $\kappa = 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,
$$\displaylines{
C'_{il} = \gamma_{il} x_i y_l - \rho_{il}
C_{il}, \cr
\hfill\llap (2.2)
\cr
D'_{il} = p_{il} \rho_{il} C_{il} - \sigma_{il} D_{il}.
}
$$
Here $1/ \gamma_{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/\rho_{il}$ is the average length of the
first stage and $1/\sigma_{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.
\begin{theorem} % 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.
\end{theorem}
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 {\it 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 {\it 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 {\it Heesterbeek} and {\it 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),
$$\displaylines{
C^*_{il} = \delta_{il} x^*_i y^*_l, \quad
D^*_{il}= {p_{il} \rho_{il} \over \sigma_{il}} C^*_{il}, \cr
\hfill x_i^* = \pi_i X_i - \sum_{k=1}^n [C^*_{ik} + \kappa D^*_{ik}],
\hfill\llap (2.3) \cr
y_l^* = \tilde \pi_l Y_l - \sum_{j=1}^m [C^*_{jl} + D^*_{jl} ],
}$$
with
$$
\delta_{il} = {\gamma_{il} \over \rho_{il}}.
\eqno(2.4)
$$
We expedite the algebra by introducing
$$\displaylines{
x_i^* = u_i \pi_i X_i, \cr
\hfill y_l^* = w_l \tilde \pi_l Y_l, \hfill\llap (2.5) \cr
C^*_{il} = K_{il} \delta_{il} \pi_i X_i \tilde \pi_l Y_l.
}$$
Then
$$\displaylines{
K_{il} = u_i w_l, \cr
\hfill u_i = 1 - \sum_{k=1}^n \tilde \pi_k Y_k \xi_{ik} \delta_{ik}
K_{ik}, \hfill \llap{(2.6)} \cr
w_l = 1 - \sum_{j=1}^m \pi_j X_j \tilde \xi_{j l} \delta_{jl} K_{jl},
}$$
with
$$
\displaylines{
\hfill \xi_{jk} = 1 + \kappa {p_{jk} \rho_{jk} \over \sigma_{jk}},
\quad \kappa \in \{0,1\},
\qquad
\tilde \xi_{jk} = 1 + {p_{jk} \rho_{jk} \over \sigma_{jk}}.
\hfill\llap(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).
\begin{theorem} % 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).
\end{theorem}
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} \ge
{ \pi_i X_i \over
1+ \sum_{j=1}^m \tilde \xi_{jl} \delta_{jl} \pi_j X_j
}
\delta_{il}
{\tilde \pi_l Y_l
\over
1 + \sum_{k=1}^n \xi_{ik} \delta_{ik} \tilde \pi_k Y_k },
\eqno(2.8)
$$
the estimate from above has the following form,
$$
C_{il}\le{\pi_i X_i\over\displaystyle{ 1+
\sum_{j=1}^m { \tilde \xi_{jl} \delta_{jl} \pi_j X_j
\over 1 +\sum_{k=1}^n \xi_{jk} \delta_{jk} \tilde \pi_k Y_k}}}
\delta_{il}{\tilde \pi_l Y_l
\over\displaystyle{ 1
+\sum_{k=1}^n { \xi_{ik} \delta_{ik} \tilde \pi_k Y_k
\over 1 + \sum_{j=1}^m \tilde \xi_{jk} \delta_{jk} \pi_j X_j}}}.
\eqno(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 $\gamma_{il}$ are very small compared with the $\rho_{il} $
and $\sigma_{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 {\it 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}
\approx
{\delta_{il} \pi_i \tilde \pi_l X_i Y_l
\over 1 + \sum_{k=1}^n \xi_{ik} \delta_{ik} \tilde \pi_k Y_k
+\sum_{j=1}^m \tilde \xi_{jl} \delta_{jl} \pi_j X_j}.
\eqno(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 $\xi_{ik} <\tilde\xi_{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} \approx
{ \pi_i X_i \over
\displaystyle{ 1+
\sum_{j=1}^m \tilde \xi_{jl} \delta_{jl} \pi_j X_j
}
}
\delta_{il}
{
\tilde \pi_l Y_l
\over
\displaystyle{ 1 +
\sum_{k=1}^n { \xi_{ik} \delta_{ik} \tilde \pi_k Y_k
\over
1 + \sum_{j=1}^m \tilde \xi_{jk} \delta_{jk} \pi_j X_j } } }.
$$
Multiplying the factors out, we obtain
$$
C_{il}
\approx
{
\pi_i X_i \delta_{il} \tilde \pi_l Y_l
\over
\displaystyle{
1+
\sum_{j=1}^m \tilde \xi_{jl} \delta_{jl} \pi_j X_j
+
\sum_{k=1}^n \xi_{ik} \delta_{ik} \tilde \pi_k Y_k
{
1+ \sum_{j=1}^m \tilde \xi_{jl} \delta_{jl} \pi_j X_j
\over
1 + \sum_{j=1}^m \tilde \xi_{jk} \delta_{jk} \pi_j X_j } } }.
\eqno(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 {\it De
Boer} and {\it Perelson} (1995) obtain for T cell proliferation by
a Pad\'e 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}
\approx
{
\delta_{il} \pi_i \tilde \pi_l X_i Y_l
\over
1 +
\max \left \{ \;\sum_{k=1}^n \xi_{ik} \delta_{ik} \tilde \pi_k Y_k
, \;
\sum_{j=1}^m \tilde \xi_{jl} \delta_{jl} \pi_j X_j
\right \}
}.
\eqno(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)),
$$
\displaylines{
\hfill \xi_{jk} \delta_{jk} = {\gamma_{jk} \over \rho_{jk} } + \kappa
{p_{jk} \gamma_{jk} \over \sigma_{jk} }, \qquad
\tilde \xi_{jk} \delta_{jk}
= {\gamma_{jk} \over \rho_{jk}} + {p_{jk} \gamma_{jk} \over \sigma_{jk} }.
\hfill\llap (2.13)}
$$
Again, under the condition that
the $\gamma_{il}$ are very small compared with the $\rho_{il} $
and $\sigma_{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).
\begin{theorem} % 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
$$\displaylines{
\theta \sum_{l=1}^n \tilde \pi_l Y_l
+\tilde \theta \sum_{i=1}^m \pi_i X_i < 1,\cr
\theta= \max_{ik}
\Big({\gamma_{ik} \over \rho_{ik} } + \kappa {p_{ik} \gamma_{ik} \over
\sigma_{ik} }\Big),\cr
\tilde \theta = \max_{ik}
\Big({\gamma_{ik} \over \rho_{ik} } + {p_{ik} \gamma_{ik} \over
\sigma_{ik} }\Big).
}$$
\end{theorem}
For local asymptotic stability alone, the assumption in Theorem 2.3
can be slightly relaxed.
\section{Proportionate mixing}
The term {\it 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
$$
{\gamma_{jl} \over \rho_{jl}} , \quad {p_{jl} \gamma_{jl} \over \sigma_{jl}}
$$
are independent of $j$ and $l$.
In other words, if there exist constants $c, \tilde c >0$ such that
$$
\gamma_{jl} = c \rho_{jl} = \tilde c {\sigma_{jl} \over 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
$ \xi_{jk} \delta_{jk} $ and $ \tilde \xi_{j k} \delta_{jk} $
do not depend on $j$ and $k$. For the subsequent considerations it is
actually sufficient that
$
\xi_{jk} \delta_{jk}
$
do not depend on $j$, while
$
\tilde \xi_{j k} \delta_{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 - \bar y K) (1 - \bar x K)
$$
with
$$
\bar x = \sum_{j=1}^m \pi_j X_j \tilde \xi_{j l} \delta_{jl},
\quad
\bar y = \sum_{k=1}^n \tilde \pi_k Y_k \xi_{ik} \delta_{ik}. \eqno(3.1)
$$
Solving this quadratic equation we obtain
$$
K={1 \over 2 \bar x \bar y}
\left
[1 + \bar x + \bar y - \sqrt {(1+ \bar x+ \bar y)^2 -
4 \bar x \bar y }
\right ].
$$
The following form is more useful,
$$
K={2\over
1 + \bar x + \bar y + \sqrt { 1 + 2 (\bar x+ \bar y) + (\bar x- \bar y)^2}}.
\eqno(3.2)
$$
Among other things it reveals why we have taken this and not the other
root which would be infinity for $\bar x = \bar y = 0$.
Recalling (2.3) we have
$$\displaylines{
\hfill D_{il}^* = {p_{il} \rho_{il} \over \sigma_{il}} C_{il}^*,
\hfill\llap(3.3)\cr
C_{il}^*= {2 \delta_{il} \pi_i X_i \tilde \pi_l Y_l
\over
1 + \bar x + \bar y + \sqrt { 1 + 2 (\bar x+ \bar y) + (\bar x- \bar y)^2}},
}$$
with $\bar x $ and $ \bar y$ from (3.1). We call this formula for
$C_{il}$ and $D_{il}$ the {\it 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. {\it 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.
\section{Estimates and comparisons}
The formulas (3.2) and (3.3) in the previous section show that the function
$$
\zeta (\bar x, \bar y) =
1 + \bar x + \bar y + \sqrt { 1 + 2 (\bar x+ \bar y) + (\bar x- \bar y)^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.
\begin{lemma} % Lemma 3.1.
$2(1 + \max \{\bar x, \bar y \} )\le
\zeta (\bar x, \bar y)\le 2 (1 + \bar x + \bar y)$.
\end{lemma}
\paragraph{Proof:} The second estimate is obvious. As for the first,
$$
\zeta(\bar x, \bar y)\ge
1 + \bar x + \bar y + \sqrt { 1 + 2 |\bar x - \bar y| + (\bar x- \bar y)^2}
=1 + \bar x + \bar y + 1 + |\bar x - \bar y|.
$$
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
$\zeta$.
\begin{lemma} % Lemma 3.2.
$$\displaylines{
\hfill {2 \bar x \bar y \over 1 + \bar x + \bar y}
\le 2(1 + \bar x + \bar y) - \zeta (\bar x, \bar y)
\le {2 \bar x \bar y \over 1 + \max \{\bar x , \bar y\}},
\hfill\llap{\rm (a)}\cr
\hfill { 2 \min \{\bar x, \bar y \} \over 1 + \max \{\bar x , \bar y\}}
\le \zeta (\bar x, \bar y ) - 2 ( 1+ \max \{\bar x, \bar y \})
\le { 2 \min \{\bar x, \bar y \} \over 1 + |\bar x - \bar y | }.
\hfill\llap{\rm (b)}
}$$
\end{lemma}
\paragraph{Proof:} Notice that
\begin{eqnarray*}
2(1 + \bar x + \bar y) - \zeta (\bar x,\bar y)
&=&1 + \bar x + \bar y -\sqrt {
1 + 2 (\bar x + \bar y) + (\bar x - \bar y )^2 } \\
&= &{4 \bar x \bar y \over1 + \bar x + \bar y + \sqrt {
1 + 2 (\bar x + \bar y) + (\bar x - \bar y )^2 }}.
\end{eqnarray*}
Replacing $\bar x + \bar y $ by $|\bar x - \bar y |$ or vice versa
yields the two estimates in (a). Further
\begin{eqnarray*}
\zeta (\bar x,\bar y) - 2 (1+ \max \{ \bar x, \bar y \} )
&=& \sqrt {1 + 2 (\bar x + \bar y) + (\bar x - \bar y )^2 }
- (1 + |\bar x- \bar y|) \\
&=& {2 (\bar x+ \bar y - |\bar x-\bar y| )
\over\sqrt {1 + 2 (\bar x + \bar y) + (\bar x - \bar y )^2 }
+ (1 + |\bar x-\bar y|) } \\
&=&{4 \min \{\bar x , \bar y \}
\over\sqrt {1 + 2 (\bar x + \bar y) + (\bar x - \bar y )^2 }
+ (1 + |\bar x-\bar y|) }
\end{eqnarray*}
Replacing $\bar x + \bar y $ by $|x - y |$ or vice versa under the
square root yields the two estimates in (b).
\subsection*{Approximations}
Lemma 3.2 (a) shows that, if $\bar x$ and $\bar y$ are small,
$\zeta (\bar x ,\bar y)\approx 2 (1 + \bar x + \bar y) $ and we can
approximate $C$ by a functional response \`a la Beddington and
DeAngelis et al.
$$
C_{il}
\approx
{\delta_{il} \pi_i X_i \tilde \pi_l Y_l
\over
1 + \bar x + \bar y }.
\eqno(3.4)
$$
The functional response (3.4) can also be derived as an
approximation of the square root formula obtained by {\it Ruxton et
al.} (1992) and from effects of spatial grouping ({\it Cosner et
al.}, 1999).
However, if one of $\bar x$ and $\bar y$ is small and the other
large, the following approximation of {\it flip-flop} type is
better than the approximation (3.4),
$$
C_{il}\approx{\delta_{il} \pi_i X_i \tilde \pi_l Y_l
\over 1 + \max \{ \bar x, \bar y \}}.\eqno(3.5)
$$
As seen from Lemma 3.2, the flip-flop response is no good
approximation at all, if $\bar x$ and $\bar y$ 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 predator-prey models, if $\bar y < \bar 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 {\it 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 {\it Heesterbeek} and {\it Metz}
(1993).
If we are dealing with one species on each side only, we have
$$
C = {2 \over \delta \xi \tilde \xi} { \bar x \bar y \over \zeta
(\bar x, \bar y ) }.
$$
If one of $\bar x$ and $\bar y$ is large and the other
small, the flip-flop approximation (3.5) can be further
approximated by
$$
C \approx {1 \over \delta \xi \tilde \xi} { \bar x \bar y \over
\max\{\bar x, \bar y \} }
= {1 \over \delta \xi \tilde \xi}
\min \{\bar x, \bar y \}
$$
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 \approx {1 \over \delta \xi \tilde \xi} { \bar x \bar y \over 1 +
\bar x + \bar y },
$$
we would need to assume that $\bar x + \bar y$ is large, while in
justifying (3.4) we needed to assume that this sum is small.
\section{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 \`a
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' = \tilde r X - \tilde \nu X^2 - p \rho C
$$
where $C$ denotes the first stage of the prey-predator complex.
Recall that $\rho C$ is the exit rate of the first stage and $p
\rho 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),
$$\displaylines{
C= {2 \delta \pi X \tilde \pi Y
\over 1 + \bar x + \bar y + \sqrt { 1 + 2 (\bar x+ \bar y)
+ (\bar x- \bar y)^2}}. \cr
\bar x = \pi X \tilde \xi \delta\,,
\quad \bar y = \tilde \pi Y \xi \delta\,.
}$$
One may wonder whether $\tilde r X$ should
be replaced by $\tilde r (X- C) $ and $\tilde \nu X^2$ by
$\tilde \nu (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 \rho C$.
Introducing new parameters and a new time
leads to the following equation for
$\bar x$:
$$
{\bar x'\over \bar x} = r - \nu \bar x -
{2 \over
1 + \bar x + \bar y + \sqrt { 1 + 2 (\bar x+ \bar y) + (\bar x- \bar y)^2}}.
\eqno(4.1)
$$
Notice that $\bar y$ 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
$$
{\bar x'\over \bar x} = r - \nu \bar x -
{1 \over
1 + \bar x + \bar y },
\eqno(4.2)
$$
and
$$
{\bar x'\over \bar x} = r - \nu \bar x -
{1 \over
1 + \max \{ \bar x , \bar y \}}.
\eqno(4.3)
$$
Equilibria of equation (4.1) are intersection points of the
decreasing line $ (1/2)(r - \nu \bar x)$ and the graph of the
function
$$
\phi (\bar x) =
{1 \over
1 + \bar x + \bar y + \sqrt { 1 + 2 (\bar x+ \bar y) + (\bar x- \bar y)^2}}.
$$
One easily sees that $\phi$ is strictly decreasing. As far as convexity
is concerned, there are two cases (Appendix E):
\begin{description}
\item{Case 1:} $\bar y \le 1$. Then $\phi $ is strictly convex.
\item{Case 2:} $\bar y > 1 $.
There exists $\bar x^\sharp > 0 $ such that $\phi$ is strictly
concave on $[0 , \bar x^\sharp] $ and
strictly convex on $[\bar x^\sharp, \infty )$.
\end{description}
\noindent
In the second case, we could not determine $\bar x^\sharp$
precisely, but at least locate it in the interval $(\bar y -1,
(5/3) \bar y - 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.
\section{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,
${\mathcal 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} = {S_i \over X_i} \sigma_{ik} q_{ik} D_{ik} {\tilde I_k
\over Y_k}. \eqno(5.1)
$$
Here $S_i$ is the number of susceptibles in female
group $i$ and ${S_i \over X_i}$ is the probability that the female
group $i$ individual involved in the complex is actually
susceptible.
$\tilde I_k$ is the number of infectives in male group $k$
and ${\tilde I_k \over 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} = {S_i \over X_i} p_{ik} \rho_{ik} q_{ik} C_{ik}
{\tilde I_k \over Y_k},
$$
and by (2.5),
$$
J_{ik}
=
\pi_i S_i p_{ik} \gamma_{ik} q_{ik} K_{ik}
\tilde \pi_k \tilde I_k,
\eqno(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
$$
\tilde J_{ik}
=
{I_i \over X_i} \sigma_{ik} q_{ik} D_{ik}
{\tilde S_k \over Y_k}
=
\pi_i I_i p_{ik} \gamma_{ik} \tilde q_{ik} K_{ik}
\tilde \pi_k \tilde S_k.
$$
The incidence in female group $i$ is now given by $\sum_{k=1}^n
J_{ik}$, while the incidence in male group $l$ is given by
$\sum_{j=1}^m \tilde 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),
$$\displaylines{
K_{il} \approx { 1 \over
1 + \sum_{k=1}^n \xi_{ik} \delta_{ik} \tilde \pi_k Y_k }
\; {1 \over 1 + \sum_{j=1}^m \tilde \xi_{jl} \delta_{jl} \pi_j X_j},
\cr
K_{il} \approx { 1 \over
1 + \sum_{k=1}^n \xi_{ik} \delta_{ik} \tilde \pi_k Y_k
+\sum_{j=1}^m \tilde \xi_{jl} \delta_{jl} \pi_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).
\subsection*{The basic reproduction ratio}
Let $\tau_i$ the average length of the infective period for
female group $i$ and $\tilde \tau_k$ the
average length of the infective period for male group $k$.
Let ${\mathcal R}_{ik} $ be the {\it 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. {\it Jacquez,
Simon, Koopman}, 1995, where ${\mathcal R}_{ik} $ is called the basic
reproduction number for the female group $i$ contacts of an
infected male in group $k$.) We have
$$
{\mathcal R}_{ik}
=
J_{ik} {\tilde \tau_k \over I_k},
$$
with $X_i$ replacing $S_i$ in (5.2), i.e.,
$$
{\mathcal R}_{ik}=
\pi_i X_i p_{ik} \gamma_{ik} q_{ik} K_{ik}
\tilde \pi_k \tilde \tau_k .
$$
Similarly we define the basic reproduction ratio from female group
$i$ to male group $k$, $\tilde {\mathcal R}_{ki}$, and obtain
$$
\tilde {\mathcal R}_{ki}
=\pi_i \tau_i p_{ik} \gamma_{ik} \tilde q_{ik} K_{ik}
\tilde \pi_k Y_k.
$$
Let ${\mathcal R}$ be the $m \times n $ matrix with entries ${\mathcal R}_{ik}$
and $\tilde {\mathcal R} $ the $n \times m $ matrix with entries $\tilde {\mathcal 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
$$
{\mathcal R}^\diamond = \pmatrix { 0 & {\mathcal R} \cr
\tilde {\mathcal R} & 0 \cr }.
$$
The basic reproduction ratio of the disease, ${\mathcal R}_0$, is mathematically
defined as the spectral radius of the basic reproduction matrix,
${\mathcal R}^\diamond$. 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).
${\mathcal R}_0$ is associated with a threshold phenomenon for the disease-free
state of the population. If ${\mathcal R}_0 < 1 $, the disease-free 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 ${\mathcal R}_0 > 1 $, the disease-free state is unstable and often the
disease becomes endemic in the population.
We have
$$
({\mathcal R}^\diamond)^2
=\pmatrix {{\mathcal R} \tilde {\mathcal R} & 0 \cr
0 & \tilde {\mathcal R} {\mathcal R} \cr }
$$
which implies that
${\mathcal R}_0^2$ is the spectral radius of the matrix
${\mathcal R} \tilde {\mathcal R} $.
(Cf. Diekmann, Dietz, Heesterbeek, 1991.)
\subsection*{The case of separable mixing}
Here we assume that the matrices $p_{ik}, \gamma_{ik}, q_{ik},
\tilde q_{ik} $ are separable, i.e., they are of the form $\alpha_i
\beta_k$. Since the matrix $K_{ik}$ is separable as well by (2.6), the
matrices ${\mathcal R}$ and $\tilde {\mathcal R}$ are also separable, i.e., of the
form
$$\displaylines{
{\mathcal R}_{jk} = \alpha_j \beta_k \cr
\tilde {\mathcal R}_{kj} = \tilde \alpha_j \tilde \beta_k\,.
}$$
The product of these two separable matrices is also separable and
its spectral radius is given by
$$
\big( \sum_{j=1}^m \alpha_j \tilde \alpha_j \Big)
\Big( \sum_{k=1}^n \beta_k \tilde \beta_k \Big).
$$
Hence
$$
{\mathcal R}_0^2 = \sum_{j=1}^m \sum_{k=1}^n {\mathcal R}_{jk} \tilde {\mathcal R}_{kj}.
$$
Substituting the expressions above,
$$
{\mathcal R}_0^2 = \sum_{j=1}^m \sum_{k=1}^n
(\pi_j)^2 (\tilde \pi_k)^2
p_{jk}^2 \gamma_{jk}^2 q_{jk} \tilde q_{jk} K_{jk}^2
\tau_j \tilde \tau_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, \ldots, X_m), Y= Y_1, \ldots, 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.
\section{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
predator-prey 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
{\bar x \bar y \over 1 + \bar x + \bar y} $, with $\bar x, \bar y$
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{\bar x \bar y \over 1 + \max
\{\bar x, \bar y\}}$
is a better approximation of the quasi-steady state
in complex formation (Section 3). If
$ \max \{\bar x, \bar y\}$ is large compared to one, the
flip-flop relationship takes the form of a familiar marriage function,
$c \min\{\bar x, \bar y \}$. 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).
\paragraph{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.
\subsection*{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 \times n $ matrices which can be identified with ${\mathbb R}^{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 \to b$ in case that $b < \infty$.
Define a closed convex set ${\cal C} $ as the set of those $C_{il} \ge 0$ such
that $x_i, y_l \ge 0$ with these being given by (2.1).
From (2.1) we obtain that, on $[0,b)$,
\begin{eqnarray*}
x_i'&=&- \sum_{k=1}^n C'_{ik} - \kappa \sum_{k=1}^n D'_{ik} \\
&= &\sum_{k=1}^n \rho_{ik} [1 - \kappa p_{ik}] C_{ik}
+\kappa \sum_{k=1}^n \sigma_{ik} D_{ik}
-\sum_{k=1}^n \gamma_{ik} x_i y_k, \\
y_l'&=&\sum_{j=1}^m \rho_{jl} [1 - p_{jl}] C_{jl}
+\sum_{j=1}^m \sigma_{jl} D_{jl}
-\sum_{j=1}^m \gamma_{jl} x_j y_l\,.
\end{eqnarray*}
It follows from Proposition B.7 in Smith, Waltman (1995)
that the convex set ${\cal C}$ is forward invariant under the solution.
By (2.1), solutions $C_{ik}\ge 0 $ with $x_i, y_l \ge 0$ are bounded,
which implies that solutions which satisfy these constraints initially
are defined for all forward times and satisfy the constraints forever.
\subsection*{Appendix B: Uniqueness of the steady state}
We will show that the system (2.6)
for given $X = (X_1, \ldots, X_m), Y = (Y_1, \ldots, 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),
$$\displaylines{
K_{il} = u_i w_l\cr
\hfill u_i = 1 - \sum_{k=1}^n \alpha_{ik} K_{ik} \hfill\llap{(B1)}\cr
w_l = 1 - \sum_{j=1}^m \tilde \alpha_{jl} K_{jl} ,
}$$
where we are interested in $u_i$ and $w_l$ being non-negative.
Here
$$
\alpha_{il} = \xi_{il} \delta_{il} \tilde \pi_l Y_l\,,
\quad
\tilde \alpha_{il} = \tilde \xi_{il} \delta_{il} \pi_i X_i\,.
\eqno(B2)
$$
In (B1) we substitute the first into the second and third equation,
$$
u_i = 1 - u_i \big( \sum_{k=1}^n \alpha_{i k} w_k \big)\,,
\quad
w_l = 1 - w_l \big( \sum_{j=1}^m \tilde \alpha_{jl} u_j \big)\,.
$$
We solve for $u_i$ and $w_l$,
$$
u_i = {1 \over 1 + \sum_{k=1}^n \alpha_{ik} w_k }\,,
\quad
w_l = {1 \over 1 + \sum_{j=1}^m \tilde \alpha_{jl} u_j }.
\eqno(B3)
$$
Notice that any non-negative solution $u = (u_1, \ldots, u_m),
w = (w_1, \ldots, w_n) $ of (B3) satisfies
$$
u_i \le 1, \quad w_l \le 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 = {1 \over 1 +
\displaystyle{\sum_{k=1}^n
{\alpha_{ik} \over 1 + \sum_{j=1}^m \tilde \alpha_{jk} u_j }
}}
\quad =: G_i (u).
\eqno(B4)
$$
Let $G : [0, \infty)^m \to [0, \infty)^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 ${\mathbb R}^m$ by setting
$$
x \le y \iff x_i \le y_i, \; i=1, \ldots, m,
$$
for vectors $ x = (x_1, \ldots, x_m), y = (y_1, \ldots, y_m )$.
We see that $G$ is an increasing map, i.e.,
$$
G (u) \le G(v) \hbox{ whenever } u \le v.
$$
We also notice that $G(u) \le {\bf 1} $ for all $u \in [0, \infty)^m$
where ${\bf 1} = (1, \ldots, 1)$.
Let $[0,{\bf 1}]$ denote the generalized order interval of those
vectors $u$ with $ 0 \le u \le {\bf 1}$. $G$ maps the closed convex bounded
set $[0, {\bf 1}]$ into itself and is continuous. Hence $G$ has a fixed
point $u$ in $[0,{\bf 1}]$ by Brouwer's fixed point theorem.
Actually $G$ maps $[0, \infty)^m$ into the order interval $ [G(0), {\bf 1}]
$,
and $ G(0) \ge \delta {\bf 1}$ for some $\delta >0$.
In order to show that there exists exactly one fixed point of $G$ in
$[0, \infty)^m$, we use the concavity (or sub-homogeneity)
technique of Krasnosel'skii (1964), Section 6.1.
\smallskip
Let $v$ be a positive vector and $\xi \in (0,1)$. Then
\begin{eqnarray*}
G_i (\xi v)
&= & {1 \over 1 + \sum_{k=1}^n \alpha_{ik}
\Big( {1 \over
1 + \sum_{j=1}^m \tilde \alpha_{jk} \xi v_j } \Big)}
\ge {1 \over 1 + \sum_{k=1}^n \alpha_{ik}
\Big( {1 \over
\xi + \sum_{j=1}^m \tilde \alpha_{jk} \xi v_j } \Big)} \\
&= & {1 \over 1 + {1 \over \xi} \sum_{k=1}^n \alpha_{ik}
\Big( {1 \over 1 + \sum_{j=1}^m \tilde \alpha_{jk} v_j } \Big)}
= {\xi \over \xi + \sum_{k=1}^n \alpha_{ik}
\Big( {1 \over
1 + \sum_{j=1}^m \tilde \alpha_{jk} v_j } \Big) }\\
&>& \xi { 1 \over 1 + \sum_{k=1}^n \alpha_{ik}
\Big( {1 \over 1 + \sum_{j=1}^m \tilde \alpha_{jk} v_j } \Big) }
= \xi G_i (v).
\end{eqnarray*}
This implies that $G$ is ${\bf 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,
$$\displaylines{
u(0) = 0, \quad u(j+1) = G (u(j)), \cr
v(0) = {\bf 1}, \quad 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(\bf 1)$ and analogous considerations for
$w$ in (B3) rather than $u$ (see also the next paragraph).
\subsection*{A Ross type estimate and approximation}
From (B3) we obtain an estimate from below for $K_{il}$,
$$
K_{il} = u_i w_l \ge
{1 \over 1 + \sum_{k=1}^n \alpha_{ik} }
\;
{1 \over
1 + \sum_{j=1}^m \tilde \alpha_{jl} }
\eqno(B5)
$$
which leads to estimate (2.8) involving an expression
reminiscent of Ross solutions (Castillo-Chavez, Busenberg, 1991).
If the $\alpha$'s and $\tilde \alpha$'s are small, (B5) actually
provides a reasonable approximation of $K_{il}$.
\subsection*{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} \approx
{1 \over 1 + \sum_{k=1}^n \alpha_{ik}
+ \sum_{j=1}^m \tilde \alpha_{jl}}. \eqno(B6)
$$
The following consideration shows that the approximation (B6) may
be even better than the Ross type approximation. By repeated use of
(B1),
\begin{eqnarray*}
{1 \over K_{il}}&=& {1 \over u_i w_l}
=\Big( 1 + \sum_{k=1}^n \alpha_{ik} w_k \Big)\;
\Big( 1 + \sum_{j=1}^m \tilde \alpha_{jl} u_j \Big)\\
& =& 1 + \sum_{k=1}^n \alpha_{ik} w_k + \sum_{j=1}^m \tilde
\alpha_{jl} u_j +
\sum_{k=1}^n \sum_{j=1}^m \alpha_{ik} w_k \tilde \alpha_{jl} u_j\\
&=& 1 + \sum_{k=1}^n \alpha_{ik}
\biggl (1 - w_k \sum_{j=1}^m \tilde \alpha_{jk} u_j\biggr )
+\sum_{j=1}^m \tilde \alpha_{jl}
\biggl(1 - u_j \sum_{k=1}^n \alpha_{jk} w_k\biggr ) \\
&&+\sum_{k=1}^n \sum_{j=1}^m \alpha_{ik} w_k\tilde \alpha_{jl} u_j.
\end{eqnarray*}
Hence
$$
K_{il}={ 1 \over
1 + \sum_{k=1}^n \alpha_{ik} + \sum_{j=1}^m \tilde \alpha_{jl}
+ \sum_{k=1}^n \sum_{j=1}^m
\left [\alpha_{ik} \tilde \alpha_{jl} -
\alpha_{ik} \tilde \alpha_{jk} -
\alpha_{jk} \tilde \alpha_{jl}\right ]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\ge {1 \over 1 + \bar x +\bar y}
\ge {1 \over 1 + \bar x} \; {1 \over 1 + \bar y}.
$$
We conclude by mentioning that,
if the $\gamma_{il}$ are very small compared with the
$\rho_{il} $
and $\sigma_{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 $\alpha_{il}$ and $\tilde \alpha_{il}$ are small
and the products $\alpha_{ik} \tilde \alpha_{jl}$ are very small.
\subsection*{A specific predator-prey approximation}
The asymmetry in a predator-prey relationship ($\kappa = 0$ in
(2.7), only the predator takes active part in the second stage of
the complex) has the consequence that $\xi_{ik}$ may be
considerably smaller than $\tilde \xi_{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 $\alpha_{il}$ is considerably smaller than $\tilde
\alpha_{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
$$\displaylines{
u_i \approx {1 \over 1 +
\sum_{k=1}^n \alpha_{ik}
\Big( {1 \over 1 + \sum_{j=1}^m \tilde \alpha_{jk} } \Big) } \cr
w_l \approx {1 \over 1 + \sum_{j=1}^m \tilde \alpha_{jl} }.
}$$
From (B1), multiplying the factors out,
$$
K_{il}\approx { 1\over
\displaystyle {1 + \sum_{j=1}^m \tilde \alpha_{jl}
+\sum_{k=1}^n \alpha_{ik} {
1 + \sum_{j=1}^m \tilde \alpha_{jl}
\over 1 + \sum_{j=1}^m \tilde \alpha_{jk}}}}.
\eqno(B7)
$$
Notice that this yields the approximation (B6) if there is only
one predator species. Going back to the original variables yields
(2.11).
\subsection*{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,
$$\displaylines{
0 = x_i + \sum_{k=1}^n [C_{ik} + \kappa D_{ik}] \cr
0 = y_l + \sum_{j=1}^m [C_{jl} + D_{jl} ] \cr
C'_{il} = \gamma_{il} x_i^* y_l
+ \gamma_{il} x_i y_l^* - \rho_{il} C_{il} \cr
D'_{il} = p_{il} \rho_{il} C_{il} - \sigma_{il} D_{il}.
}$$
Substituting the first two equations into the last two,
we obtain the following linear ODE system,
$$\displaylines{
C'_{il} = - \gamma_{il} x_i^* \sum_{j=1}^m [C_{jl} + D_{jl} ]
- \gamma_{il} y_l^* \sum_{k=1}^n [C_{ik} + \kappa D_{ik}]
- \rho_{il} C_{il} \cr
D'_{il} = p_{il} \rho_{il} C_{il} - \sigma_{il} D_{il}.
}$$
The algebraic system for the associated eigenvalues $\lambda$ and
their eigenvectors has the form
$$\displaylines{
\lambda C_{il} = - \gamma_{il} x_i^* \sum_{j=1}^m [C_{jl} + D_{jl} ]
- \gamma_{il} y_l^*\sum_{k=1}^n [C_{ik} + \kappa D_{ik}]
- \rho_{il} C_{il} \cr
\lambda D_{il} = p_{il} \rho_{il} C_{il} - \sigma_{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
$$
{\lambda + \rho_{il} \over \gamma_{il} } C_{il}
=-x_i^* \sum_{j=1}^m C_{jl}
\left [1 + {p_{jl} \rho_{jl} \over \lambda + \sigma_{jl}} \right ]
-y_l^*\sum_{k=1}^n C_{ik}
\left [ 1 + \kappa {p_{ik} \rho_{ik} \over \lambda + \sigma_{ik}} \right ]
$$
where at least one $C_{il}$ is different from 0.
Taking absolute values,
$$
{|\lambda + \rho_{il}| \over \gamma_{il}} |C_{il} |
\le
x_i^* \sum_{j=1}^m |C_{jl} |
\left [1 + {p_{jl} \rho_{jl} \over |\lambda + \sigma_{jl} |} \right ]
+y_l^*\sum_{k=1}^n | C_{ik} |
\left [ 1 + \kappa {p_{ik} \rho_{ik} \over |\lambda + \sigma_{ik}| } \right
].
$$
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 $\Re \lambda \ge 0$. Then
$$
{\rho_{il} \over \gamma_{il}} | C_{il} |
\le x_i^* \sum_{j=1}^m | C_{jl} |
\left [1 + {p_{jl} \rho_{jl} \over \sigma_{jl}} \right ]
+y_l^*\sum_{k=1}^n | C_{ik} |
\left [ 1 + \kappa {p_{ik} \rho_{ik} \over \sigma_{ik} } \right ].
$$
Set $\xi = \max_{il} {|C_{il} | \over C^*_{il} }$. Then $\xi > 0$.
Choose $i,l$ such that $ |C_{il} | = \xi C^*_{il} $. For such a pair
$i,l$ we have
$$
{1 \over \delta_{il}} \xi C^*_{il}
\le x_i^* \sum_{j=1}^m \xi C^*_{jl}
\left [1 + {p_{jl} \rho_{jl} \over \sigma_{jl}} \right ]
+ y_l^*\sum_{k=1}^n \xi C^*_{ik}
\left [ 1 + \kappa {p_{ik} \rho_{ik} \over \sigma_{ik} } \right ].
$$
Dividing by $\xi $ and using that $C^*_{il} = \delta_{il} x^*_i y^*_l$,
$$
x^*_i y^*_l
\le x_i^* \sum_{j=1}^m \delta_{jl} x^*_j y^*_l
\left [1 + {p_{jl} \rho_{jl} \over \sigma_{jl}} \right ]
+y_l^*\sum_{k=1}^n \delta_{ik} x^*_i y^*_k
\left [ 1 + \kappa {p_{ik} \rho_{ik} \over \sigma_{ik} } \right ].
$$
We divide by $x^*_i$ and $y^*_l$ and use (2.7),
$$
1\le \sum_{j=1}^m \delta_{jl} x^*_j
\tilde \xi_{jl}
+\sum_{k=1}^n \delta_{ik} y^*_k \xi_{ik} .
$$
Using transformation (2.5),
$$
1\le \sum_{j=1}^m \delta_{jl} \pi_j X_j u_j \tilde \xi_{jl}
+\sum_{k=1}^n \delta_{ik} \tilde \pi_k Y_k w_k\xi_{ik}.
$$
Since $u_j, w_k < 1 $, we obtain a contradiction by assuming
$$
\sum_{j=1}^m \delta_{jl} \pi_j X_j
\tilde \xi_{jl}
+\sum_{k=1}^n \delta_{ik} \tilde \pi_k Y_k
\xi_{ik}\le 1
\quad \forall i= 1, \ldots, m, \; l = 1, \ldots, n.
\eqno(C1)
$$
So this is a sufficient condition for local asymptotic stability of
the steady state.
Retracing the definition of $\alpha_{il}$ and $\tilde \alpha_{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 \le {1 \over 1 +
\sum_{k=1}^n \alpha_{ik}
\Big( {1 \over
1 + \sum_{j=1}^m \tilde \alpha_{jk} } \Big)
}
$$
and a similar estimate for $w_l$.
\subsection*{Appendix D: Global asymptotic stability of the steady state}
If $f:[0, \infty) \to {\mathbb R} $ is bounded, let
$$
f^\infty = \limsup_{t \to \infty} f(t),
\quad
f_\infty = \liminf_{t \to \infty} f(t).
$$
By the Fluctuation Lemma ({\it Hirsch et al.}, 1995)
there exists a sequence $t_j \to \infty$ (which depends
on $i$ and $l$) such that
$$
C_{il} (t_j) \to C_{il}^\infty, \quad
C_{il}'(t_j) \to 0, \quad j \to \infty.
$$
The first equation in (2.2) implies the estimate
$$
C_{il}^\infty \le \delta_{il} x_i^\infty y_l^\infty,
\eqno(D1)
$$
with $\delta_{ij}$ in (2.4).
Applying the same procedure to the second equation in (2.2),
$$
D_{il}^\infty \le {p_{il} \rho_{il} \over \sigma_{il}} C_{il}^\infty.
\eqno(D2)
$$
Using analogous arguments for the limit inferior we have
$$
C_{il \infty} \ge \delta_{il} x_{i \infty } y_{l \infty}\,,
\quad
D_{il \infty } \ge {p_{il} \rho_{il} \over \sigma_{il}} C_{il \infty}.
\eqno(D3)$$
Subtracting the first inequality in (D3) from (D1),
$$
0 \le C_{il}^\infty - C_{il \infty}
\le
\delta_{il} x_i^\infty (y_l^\infty - y_{l \infty} )
+ \delta_{il} (x_i^\infty - x_{i \infty}) y_{l \infty}.
$$
From (2.1),
$$
0
\le
C_{il}^\infty - C_{il \infty}
\le\delta_{il} \pi_i X_i (y_l^\infty - y_{l \infty} )
+\delta_{il} \tilde \pi_l Y_l (x_i^\infty - x_{i \infty}).
\eqno(D4)
$$
Further, from (2.1),
$$\displaylines{
x_i^\infty\le \pi_i X_i -
\sum_{k=1}^n [C_{ik \infty} + \kappa D_{ik \infty} ]
\le \pi_i X_i - \sum_{k=1}^n \xi_{ik} C_{ik \infty} ,
\cr
y_l^\infty \le \tilde \pi_l Y_l
- \sum_{j=1}^m \tilde \xi_{jl} C_{jl \infty} ,
\quad
x_{i \infty } \ge \pi_i X_i - \sum_{k=1}^n \xi_{ik} C_{ik}^\infty,
\quad
y_{l \infty } \ge \tilde \pi_l Y_l
- \sum_{j=1}^m \tilde \xi_{jl} C_{jl}^\infty.
}$$
Substituting these inequalities into (D4),
$$
0\le C_{il}^\infty - C_{il \infty}
\le \delta_{il} \pi_i X_i
\sum_{j=1}^m \tilde \xi_{jl} [ C_{jl}^\infty - C_{jl \infty} ]
+\delta_{il} \tilde \pi_l Y_l
\sum_{k=1}^n \xi_{ik} [ C_{ik}^\infty - C_{ik \infty} ].
$$
Set $\displaystyle C_{il}^\infty - C_{il \infty}=
\delta_{il} \Delta_{il}$,
then
$$
0\le\Delta_{il}\le\tilde \pi_l Y_l
\sum_{k=1}^n \delta_{ik} \xi_{ik} \Delta_{ik}
+\pi_i X_i \sum_{j=1}^m \delta_{jl} \tilde \xi_{jl} \Delta_{jl}.
$$
Set $\displaystyle \bar \Delta = \sum_{i=1}^m \sum_{l=1}^n \Delta_{il}$,
then
$$
0 \le \bar \Delta
\le
\left (
\theta \sum_{l=1}^n \tilde \pi_l Y_l
+\tilde \theta \sum_{i=1}^m \pi_i X_i
\right )
\bar \Delta,
\eqno(D5)
$$
with
$$\displaylines{
\theta= \max \; \left \{ \delta_{ik} \xi_{ik};
\;i =1, \ldots, m, k =1, \ldots, n \right \}, \cr
\tilde \theta = \max \; \left \{ \delta_{ik} \tilde \xi_{ik};
\;i =1, \ldots, m, k =1, \ldots, n \right \}.
}$$
By (D5), if
$$
\theta\sum_{l=1}^n \tilde \pi_l Y_l
+\tilde \theta\sum_{i=1}^m \pi_i X_i < 1,\eqno(D6)
$$
$\bar \Delta = 0 $. Since all the summands in $\bar \Delta $ are
non-negative, $ \Delta_{ik} = 0 $ for all $i,k$ and so $
C_{ik}^\infty - C_{ik \infty} = 0 $. This implies that $C_{ik} (t)
$ converges as $t \to \infty$, 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.
\subsection*{Appendix E: Proof of cases 1 and 2 in Section 4}
Recall
$$
\displaylines{
\phi (\bar x) = {1 \over
1 + \bar x + \bar y + \sqrt {\psi (\bar x )}},
\cr
\hfill\llap{(E1)}
\cr
\psi(\bar x) = 1 + 2 (\bar x+ \bar y) + (\bar x- \bar y)^2
=( 1 + \bar x - \bar y)^2 + 4 \bar y\,.
}$$
Taking the first and second derivative,
$$
\phi' (\bar x)
=-{ 1 + (\psi (\bar x ))^{-1/2} (1 + \bar x - \bar y)
\over\Big(1 + \bar x + \bar y + (\psi (\bar x ))^{1/2} \Big)^2},
$$
\begin{eqnarray*}
- \phi'' (\bar x)
&\simeq & \left [
-(\psi(\bar x))^{-3/2} (1 + \bar x - \bar y)^2
+ (\psi (\bar x ))^{-1/2}\right ]
\Big(1 + \bar x + \bar y + (\psi (\bar x ))^{1/2} \Big) \\
&& -2 [1 + (\psi (\bar x ))^{-1/2} (1 + \bar x - \bar y)]^2.
\end{eqnarray*}
Here $\simeq$ means that the two expressions it connects have the
same sign. Multiplying by $(\psi(\bar x))^{3/2}$,
\begin{eqnarray*}
- \phi'' (\bar x)
&\simeq &\left [-(1 + \bar x - \bar y)^2 + \psi (\bar x) \right ]
\Big(1 + \bar x + \bar y + (\psi (\bar x ))^{1/2} \Big)\\
&& -2 (\psi(\bar x))^{1/2}
\left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y \right ]^2.
\end{eqnarray*}
Since
$- (1 + \bar x - \bar y)^2 + \psi(\bar x) = 4 \bar y$
by (E1), we have
\begin{eqnarray*}
- \phi'' (\bar x)
&\simeq & 2 \bar y
\Big(1 + \bar x + \bar y + (\psi (\bar x ))^{1/2} \Big) \\
&& - (\psi(\bar x))^{1/2}
\left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y \right ]^2
=: \chi (\bar x).
\end{eqnarray*}
It is easy to see that
$\chi (\bar x) \to - \infty $ as $\bar x \to \infty$, i.e.,
$\phi$ is convex for large $\bar x$. Further
$$
\displaylines{
\chi(0)= 4 \bar y (1 + \bar y ) -4 (1 + \bar y)
=4 (1 + \bar y) (\bar y -1 ),
\cr
\hfill\llap{(E2)}
\cr
\chi (\bar y -1)=2 \bar y ( 2 \bar y + \sqrt{4 \bar y})
- \sqrt{4 \bar y} 4 \bar y= 4 \bar y ( \bar y - \sqrt{ \bar y} ),
}$$
so $\chi (0)$ and $\chi(\bar y -1)$ have the same sign as $\bar y -1$.
So, if $\bar y >1$, $\phi$ is concave in neighborhoods of 0
and $\bar y - 1$.
\paragraph{Lemma.} % E.1
{\sl $\phi''$ has no zero in $[2 \bar y -1,\infty)$ and
$\phi$ is strictly convex on $[2\bar y -1, \infty)$. }
\paragraph{Proof:} $\phi''$ has a zero $\bar x$ if and
only if $\chi$ has the same zero $\bar x$, i.e., if and only if
\begin{eqnarray*}
\lefteqn{ 2 \bar y \Big(1 + \bar x + \bar y +
(\psi (\bar x ))^{1/2} \Big) }\\
&= & (\psi(\bar x))^{1/2}
\left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y \right ]^2 \\
&= &(\psi(\bar x))^{1/2}
\left [ (\psi (\bar x )) + 2 (\psi(\bar x))^{1/2}(1 + \bar x - \bar y)
+(1 + \bar x - \bar y)^2 \right ] \\
&= &(\psi(\bar x))^{1/2}
\Bigl ( 2(1 + \bar x - \bar y)^2 + 4 \bar y \Bigr )
+2 \Bigl ( (1+ \bar x- \bar y)^2 + 4\bar y \Bigr ) (1 + \bar x - \bar y).
\end{eqnarray*}
Here we have used (E1) several times. Reorganizing terms we find
that $\phi'' (\bar x) = 0$ if and only if
$$
(\psi(\bar x))^{1/2}
\Bigl ( 2 (1 + \bar x - \bar y)^2 + 2 \bar y \Bigr )
= 2 \bar y (1 + \bar x + \bar y )
-2 \Bigl ( (1+ \bar x- \bar y)^2 + 4 \bar y\Bigr ) (1 + \bar x - \bar y).
$$
Simplifying
\begin{eqnarray*}
\lefteqn{ (\psi(\bar x))^{1/2}
\Bigl ( (1 + \bar x - \bar y)^2 + \bar y \Bigr ) }\\
&=&- (1+ \bar x- \bar y)^3
-\bar y \Bigl ( 4 (1 + \bar x - \bar y) - (1 + \bar x + \bar y ) \Bigr )
\\
&=&-(1+ \bar x- \bar y)^3-\bar y \bigl ( 3 (\bar x + 1)
- 5 \bar y \bigr ).
\end{eqnarray*}
The left hand side of this equation is positive. Since the
right hand side is negative for $ \bar x \ge (5/3) \bar y -1 $,
there is no equality for $\bar x \ge (5/3) \bar y -1 $ and hence no
zero of $\phi''$ with $\bar x \ge (5/3) \bar y -1 $.
Since $\phi''$ is strictly positive for large $\bar
x$, we conclude that $\phi''$ is strictly positive
on $[(5/3) \bar y - 1, \infty) $.
\smallskip
Since $\phi''(\bar y -1 ) < 0$ by (E.2), $\phi''$
must change sign somewhere between $\bar y -1 $ and $(5/3) \bar y
-1$. It follows from the next lemma
that this is the only sign change of $\phi''$.
This finishes the proof of case 2 in Section 4.
\paragraph{Lemma.} % E.2
{\sl $\chi$ is strictly decreasing on $[0, \infty)$. }
\paragraph{Proof:} We first show that $\chi$ is strictly decreasing
on $[\bar y -1, \infty)$. To this end we consider
\begin{eqnarray*}
\chi'(\bar x)
&=& 2 \bar y
\Big(1 + (\psi (\bar x ))^{-1/2} (1+ \bar x- \bar y) \Big) \\
&& -(\psi(\bar x))^{-1/2} (1 + \bar x - \bar y)
\left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y \right ]^2 \\
&& -(\psi(\bar x))^{1/2}
2 \left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y \right ]
\Big( (\psi(\bar x))^{-1/2} (1 + \bar x - \bar y) + 1 \Big).
\end{eqnarray*}
Multiplying by $(\psi(\bar x))^{1/2}$ and reorganizing the terms,
\begin{eqnarray*}
\chi'(\bar x)&\simeq &
2 \bar y \Big((\psi (\bar x ))^{1/2} + 1+ \bar x- \bar y \Big)
-\left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y \right ]^3\\
&&- (\psi(\bar x))^{1/2}
\left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y) \right ]^2.
\end{eqnarray*}
Since
$$\displaylines{
(\psi (\bar x ))^{1/2} \ge 1 + |\bar x - \bar y|, \cr
\psi (\bar x ))^{1/2} + 1+ \bar x- \bar y \ge
2 + |\bar x- \bar y| + \bar x - \bar y
= 2 ( 1 + [\bar x - \bar y]_+ ),
}$$
we have
$$
\chi'(\bar x) \simeq 2 \bar y -
\left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y \right ]^2
- (\psi(\bar x))^{1/2}
\left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y) \right ].
\eqno(E.3)
$$
We see that $\chi'(\bar x) < 0$ if $ \bar x \ge \bar y -1 $.
If $\bar y \le 1 $, $\chi(0) \le 0$ by (E2) and so
that $\chi(\bar x) < 0 $ for all $\bar x > 0$.
This is case 1 in Section 4.
By (E3), $\chi'$ has a zero $\bar x$ if and only if
\begin{eqnarray*}
2 \bar y &= &
\left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y \right ]^2
+ (\psi(\bar x))^{1/2}
\left [ (\psi (\bar x ))^{1/2} + 1 + \bar x - \bar y) \right ]\\
&= &
2 \psi (\bar x) + 3 (\psi (\bar x ))^{1/2} (1 + \bar x - \bar y)
+ (1 + \bar x - \bar y)^2 \\
&= &3 (\psi (\bar x ))^{1/2} (1 + \bar x - \bar y)
+2 \bigl ( (1+ \bar x - \bar y)^2 + 4 \bar y \bigr )
+ (1 + \bar x - \bar y)^2 \\
&= &
3 (\psi (\bar x ))^{1/2} (1 + \bar x - \bar y)
+3 (1 + \bar x - \bar y)^2
+8 \bar y.
\end{eqnarray*}
Reorganizing terms, we see that $\chi'$ has a zero $\bar x$ if and only if
$$
(\psi (\bar x ))^{1/2} (1 + \bar x - \bar y)
=- (1 + \bar x - \bar y)^2-2 \bar y.
$$
Squaring both sides we obtain
$$
\psi(\bar x) (1 + \bar x - \bar y)^2
= (1 + \bar x - \bar y)^4+4 \bar y (1 + \bar x - \bar y)^2
+4 \bar y^2.\eqno(E.4)
$$
Remembering (E.1), $\psi(\bar x) =(1 + \bar x - \bar y)^2 + 4 \bar y$,
the equality (E4) can only hold if $\bar y = 0$. So $\chi'$
has no zero for $\bar y >0$. Since $\chi'$ is strictly
negative for large $\bar x$, it is strictly negative
everywhere, if $\bar y >0$.
\begin{thebibliography}{00}
\bibitem{a1} Arditi, R.; Ginzburg, L.R. (1989): Coupling in predator-prey
dynamics: Ratio dependence. {\sl J. Theor. Biol.} {\bf 139},
311-326
\bibitem{b1} Beddington, J.R. (1975): Mutual interference between parasites
or predators and its effects on searching efficiency.
{\sl J. Anim. Ecol.} {\bf 51}, 597-624
\bibitem{b2} Borghans, J.A.M.; De Boer, R.J.; Segel, L.A. (1996):
Extending the quasi-steady state approximation by changing
variables. {\sl Bull. Math. Biol.} {\bf 58}, 43-63
\bibitem{c1} Castillo-Chavez, C.; Busenberg, S. (1991): On the solution
of the two-sex mixing problem. {\sl Differential Equations
Models in Biology, Epidemiology and Ecology} (S. Busenberg,
M. Martelli; eds.), 80-98. Lecture Notes in Biomathematics {\bf 92}.
Springer
\bibitem{c2} Castillo-Chavez, C.; Velasco-Hernandez, J.X.; Fridman, S.
(1994): Modeling contact structures in biology. {\sl Frontiers
in Mathematical Biology} (S.A. Levin, ed.), 454-491. Lecture
Notes in Biomathematics {\bf 100}, Springer
\bibitem{c3} Cosner, C.; DeAngelis, D.L.; Ault, J.S.; Olson, D.B.
(1999): Effects of spatial grouping on the functional response of
predators. {\sl Theor. Pop. Biol.} {\bf 56}, 65-75
\bibitem{d1} DeAngelis, D.L.; Goldstein, R.A.; O'Neill, R. V. (1975): A model
for trophic interaction. {\sl Ecology} {\bf 56}, 881-892
\bibitem{d2} De Boer, R.J.; Perelson, A.S. (1995): Towards a general
function describing T cell proliferation.
{\sl J. theor. Biol.} {\bf 175}, 567-576
\bibitem{d3} Diekmann, O.; Dietz, K.; Heesterbeek, J.A.P. (1991):
The basic reproduction ratio for sexually transmitted diseases:
I. Theoretical considerations. {\sl Math. Biosci.} {\bf 107}, 325-339
\bibitem{d4} Diekmann, O.; Heesterbeek, J.A.P.; Metz, J.A.J. (1990):
On the definition and the computation of the basic reproduction
ratio ${\mathcal R}_0$ in models for infectious diseases in heterogeneous
populations. {\sl J. Math. Biol.} {\bf 28}, 365-382
\bibitem{d5} Diekmann, O.; Heesterbeek, J.A.P.; Metz, J.A.J. (1995):
The legacy of Kermack and McKendrick.
{\sl Epidemic Models: Their Structure and Relation to Data}
(D. Mollison, ed.), 95-115. Cambridge Univ. Press
\bibitem{e1} Edelstein-Keshet, L. (1987): {\sl Mathematical Models in Biology}.
McGraw-Hill
\bibitem{h1} Hadeler, K.P.; Waldst\"atter, R.; W\"orz-Busekros, A. (1988):
Models for pair formation in bisexual populations.
{\sl J. Math. Biol.} {\bf 26}, 635-649
\bibitem{h2} Hale, J.K. (1980): {\sl Ordinary Differential Equations}.
Krieger
\bibitem{h3} Heesterbeek, J.A.P.; Metz, J.A.J. (1993): The saturating contact
rate in marriage and epidemic models. {\sl J. Math. Biol.} {\bf 31},
529-539
\bibitem{h4} Hethcote, H.W.; Thieme, H.R. (1985):
Stability of the endemic equilibrium in epidemic models with
subpopulations. {\sl Math. Biosciences} {\bf 75}, 205-227
\bibitem{h5} Hirsch, W.M.; Hanisch, H.; Gabriel, J.P. (1985):
Differential equations models for some parasitic infections;
methods for the study of asymptotic behavior.
{\sl Comm. Pure Appl. Math.} {\bf 38}, 733-753
\bibitem{h6} Holling, C.S. (1965): The functional response of predators to prey
density and its role in mimicry and population regulation.
{\sl Mem. Ent. Soc. Canada} {\bf 45}, 1-60
\bibitem{h7} Holling, C.S. (1966): The functional response of invertebrate predators
to prey density. {\sl Mem. Ent. Soc. Canada} {\bf 48}
\bibitem{h8} Hoppensteadt, F. (1974): Asymptotic stability in singular
perturbation problems. II: Problems having matched asymptotic expansion
solutions. {\sl J. Diff. Eqn.} {\bf 15}, 510-521
\bibitem{h9} Hsu, P.-H.; Fredrickson, A.G. (1975): Population-changing
processes and the dynamics of sexual populations.
{\sl Math. Biosci.} {\bf 26}, 55-78
\bibitem{h10} Huisman, G.; De Boer, R.J. (1997): A formal derivation of the
`Beddington' functional response.
{\sl J. theor. Biol.} {\bf 185}, 389-400
\bibitem{h11} Jacquez, J.A.; Simon, C.P.; Koopman, J. (1995):
Core groups and the ${\mathcal R}_0$s for subgroups in heterogeneous SIS and SI
models. {\sl Epidemic Models: Their Structure and Relation to Data}
(D. Mollison, ed.), 279-301. Cambridge Univ. Press
\bibitem{k1} Krasnosel'skii, M.A. (1964): {\sl Positive Solutions of Operator
Equations}. Noordhoff
\bibitem{k2} Kuang, Y.; Beretta, E. (1998): Global qualitative analysis
of a ratio-dependent predator-prey system.
{\sl J. Math. Biol.} {\bf 36}, 389-406
\bibitem{l1} Lin, C.C.; Segel, L.A. (1974, 1988):
{\sl Mathematics Applied to Deterministic Problems in the Natural
Sciences}. Siam
\bibitem{m1} Michaelis, L.; Menten, M.I. (1913): Die Kinetik der
Invertinwirkung. {\sl Biochem. Z.} {\bf 49}, 333-369
\bibitem{m2} Monod, J. (1942): {\sl Recherches sur la croissance des
cultures bacteriennes.} Herman
\bibitem{m3} Monod, J. (1950): La technique de culture continue, th\'eorie
et applications. {\sl Annales de l'Institute Pasteur} {\bf 79},
390-401
\bibitem{r1} Ruxton, G.D.; Gurney, W.S.C; de Roos, A.M. (1992):
Interference and generation cycles.
{\sl Theor. Pop. Biol.} {\bf 42}, 235-253
\bibitem{s1} Smith, H.L.; Waltman, P. (1995): {\sl The Theory of Chemostat.
Dynamics of Microbial Competition}. Cambridge Univ. Press
\bibitem{t1} Thieme, H.R. (1985):
Local stability in epidemic models for heterogeneous
populations. {\sl Mathematics in Biology and Medicine} (V. Capasso,
E. Grosso, S.L. Paveri-Fontana; eds.), 185-189.
Lecture Notes in Biomathematics {\bf 57}. Springer
\end{thebibliography}
\noindent{\sc Horst R. Thieme \& Jinling Yang}\\
Department of Mathematics,
Arizona State University \\
Tempe, AZ 85287-1804, USA \\
e-mail: h.thieme@asu.edu
\end{document}