\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}