A Nutrient-Prey-Predator Model: Stability and Bifurcations

TL;DR

This paper analyzes a nutrient-prey-predator model in a chemostat, examining how nutrient levels influence system stability and bifurcations, including Hopf bifurcations and limit cycles, with analytical and numerical methods.

Contribution

It provides a comprehensive analysis of stability and bifurcations in a nutrient-prey-predator system with general functional responses, including analytical verification and numerical illustrations.

Findings

Identification of conditions for Hopf bifurcation

Existence of limit cycles under certain nutrient levels

Analytical and numerical validation of bifurcation scenarios

Abstract

In this paper we consider a model of a nutrient-prey-predator system in a chemostat with general functional responses, using the input concentration of nutrient as the bifurcation parameter. We study the changes in the existence of isolated equilibria and in their stability, as well as the global dynamics, as the nutrient concentration varies. The bifurcations of the system are analytically verified and we identify conditions under which an equilibrium undergoes a Hopf bifurcation and a limit cycle appears. Numerical simulations for specific functional responses illustrate the general results.

Figures

Equations

$$dN/dt$$

$$dP/dt$$

$$dZ/dt$$

$$f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}=D_{1}/\gamma_{1}\mbox{ and }f_{2}\bigl{(}\lambda_{Z}(D_{2})\bigr{)}=D_{2}/\gamma_{2}.$$

$$\lambda_{P}^{\prime}(D_{1})=\bigl{(}\gamma_{1}\cdot f_{1}^{\prime}(\lambda_{P}(D_{1}))\bigr{)}^{-1}\quad\text{and}\quad\lambda_{Z}^{\prime}(D_{2})=\bigl{(}\gamma_{2}\cdot f_{2}^{\prime}(\lambda_{P}(D_{2}))\bigr{)}^{-1}.\qed$$

$$U^{\prime }(t)=D\mu -DN(t)-{\gamma }_1^{-1}{D}_{1}P(t)-({\gamma }_1{\gamma }_2{)}^{-1}{D}_{2}Z(t)\le D\mu -\hat{D}U(t),$$

$$U(t)$$

$$\le \{\begin{array}{ll}U(0),& \text{if U(0)>(Dμ)/D,}\\ (D\mu )/\hat{D},& \text{if U(0)≤(Dμ)/D.}\end{array}$$

$$P(\mu,D_{1})=\bigl{(}\mu-\lambda_{P}(D_{1})\bigr{)}\frac{D\gamma_{1}}{D_{1}}.$$

$$E_{2}(\mu,D_{1},D_{2})=\bigl{(}N(\mu,D_{1},D_{2}),\lambda_{Z}(D_{2}),Z(\mu,D_{1},D_{2})\bigr{)}$$

$${\mu }_{c_1}(D_1,D_2)={\lambda }_P(D_1)+\frac{D_1{\lambda }_Z(D_2)}{D{\gamma }_1}$$

$$\displaystyle\bigl{(}\mu-N(\mu,D_{1},D_{2})\bigr{)}D-\lambda_{Z}(D_{2})f_{1}\bigl{(}N(\mu,D_{1},D_{2})\bigr{)}=0,$$

$$\displaystyle\gamma_{1}\lambda_{Z}(D_{2})f_{1}\bigl{(}N(\mu,D_{1},D_{2})\bigr{)}-D_{1}\lambda_{Z}(D_{2})-Z(\mu,D_{1},D_{2})f_{2}\bigl{(}\lambda_{Z}(D_{2})\bigr{)}=0.$$

$$0={\gamma }_{1}P{f}_1(N)-D_{1}P,$$

$$0=\bigl{(}\mu-\lambda_{P}(D_{1})\bigr{)}D-P\frac{D_{1}}{\gamma_{1}},\quad\text{with solution}\quad P=\bigl{(}\mu-\lambda_{P}(D_{1})\bigr{)}\frac{D\gamma_{1}}{D_{1}}\mathrel{\mathop{:}}=P(\mu,D_{1}).$$

$$Z(\mu,D_{1},D_{2})=\frac{\gamma_{2}}{D_{2}}\lambda_{Z}(D_{2})\bigl{(}\gamma_{1}f_{1}\bigl{(}N(\mu,D_{1},D_{2})\bigr{)}-D_{1}\bigr{)}.$$

$$0=\bigl{(}\mu-\lambda_{P}(D_{1})\bigr{)}D-\lambda_{Z}(D_{2})\frac{D_{1}}{\gamma_{1}}\;\text{with solution}\;\mu=\lambda_{P}(D_{1})+\frac{D_{1}\lambda_{Z}(D_{2})}{D\gamma_{1}}\mathrel{\mathop{:}}=\mu_{c_{1}}(D_{1},D_{2}).$$

$$J=\left[\begin{array}{ccc}-D& -f_1(N)& 0\\ {\gamma }_{1}P{f}_1^{\prime }(N)& {\gamma }_1{f}_1(N)-D_1-Z{f}_2^{\prime }(P)& -f_2(P)\\ 0& {\gamma }_{2}Z{f}_2^{\prime }(P)& {\gamma }_2{f}_2(P)-D_2\end{array}\right].$$

$$J(E_0)=\left[\begin{array}{ccc}-D& -f_1(\mu )& 0\\ 0& {\gamma }_1{f}_1(\mu )-D_1& 0\\ 0& 0& -D_2\end{array}\right],$$

$$L(N,Z)=\frac{(\mu -N{)}^2}{2}+\frac{Z^2}{2}.$$

$$L^{\prime }(N,Z)=-D(\mu -N{)}^2-D_2{Z}^2<0$$

$$J\bigl{(}E_{1}(\mu,D_{1})\bigr{)}=\\ \begin{bmatrix}-D-P(\mu,D_{1})f_{1}^{\prime}\bigl{(}\lambda_{P}(D_{1})\bigr{)}&-f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}&0\\ \gamma_{1}P(\mu,D_{1})f_{1}^{\prime}\bigl{(}\lambda_{P}(D_{1})\bigr{)}&0&-f_{2}\bigl{(}P(\mu,D_{1})\bigr{)})\\ 0&0&\gamma_{2}f_{2}\bigl{(}P(\mu,D_{1})\bigr{)}-D_{2}\end{bmatrix}.$$

$$L(N,P)=\\ \int_{\lambda_{P}(D_{1})}^{N}\frac{f_{1}(n)-f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}}{f_{1}(n)}\,dn+\frac{1}{\gamma_{1}}\Bigl{(}P-P(\mu,D_{1})-P(\mu,D_{1})\ln\bigl{(}P/P(\mu,D_{1})\bigr{)}\Bigr{)}.$$

$$\frac{\partial L}{\partial N}=\frac{f_{1}(N)-f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}}{f_{1}(N)}=0\quad\text{and}\quad\frac{\partial L}{\partial P}=\frac{1}{\gamma_{1}}\bigl{(}1-P(\mu,D_{1})/P\bigr{)}=0$$

$$\displaystyle\frac{\partial^{2}L}{\partial N^{2}}\bigl{(}\lambda_{P}(D_{1}),P(\mu,D_{1})\bigr{)}$$

$$\displaystyle\frac{\partial^{2}L}{\partial P^{2}}\bigl{(}\lambda_{P}(D_{1}),P(\mu,D_{1})\bigr{)}$$

$$f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}=\frac{D_{1}}{\gamma_{1}}\quad\text{and}\quad P(\mu,D_{1})=\frac{\bigl{(}\mu-\lambda_{P}(D_{1})\bigr{)}D}{f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}},$$

$$\begin{split}L^{\prime}(N,P)&=\frac{f_{1}(N)-f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}}{f_{1}(N)}\bigl{(}(\mu-N)D-Pf_{1}(N)\bigr{)}\\ &\quad\quad\quad\quad+\frac{1}{\gamma_{1}}\biggl{(}1-\frac{P(\mu,D_{1})}{P}\biggr{)}(\gamma_{1}f_{1}(N)-D_{1})P\\ &=\bigl{(}f_{1}(N)-f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}\bigr{)}\biggl{(}\frac{(\mu-N)D}{f_{1}(N)}-\frac{\bigl{(}\mu-\lambda_{P}(D_{1})\bigr{)}D}{f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}}\biggr{)}\\ &=\bigl{(}f_{1}(N)-f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}\bigr{)}\biggl{(}\frac{(\mu-N)}{\bigl{(}\mu-\lambda_{P}(D_{1})\bigr{)}}-\frac{f_{1}(N)}{f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}}\biggr{)}\frac{\bigl{(}\mu-\lambda_{P}(D_{1})\bigr{)}D}{f_{1}(N)}.\end{split}$$

$$\frac{\mu{-}N}{\mu{-}\lambda_{P}(D_{1})}>1>\frac{f_{1}(N)}{f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}},\;\text{and}\;\frac{\mu{-}N}{\mu{-}\lambda_{P}(D_{1})}-\frac{f_{1}(N)}{f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}}>0.$$

$$\frac{\mu-N}{\mu-\lambda_{P}(D_{1})}<1<\frac{f_{1}(N)}{f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}},$$

$$E_2({\mu }_{c_1}(D_1,D_2),D_1,D_2)=({\lambda }_P(D_1),{\lambda }_Z(D_2),0)=E_1({\mu }_{c_1}(D_1,D_2),D_1).$$

$$J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)}=\\ \begin{bmatrix}-D{-}\lambda_{Z}(D_{2})f_{1}^{\prime}(N(\mu,\!D_{1},\!D_{2}))&\mspace{-15.0mu}-f_{1}(N(\mu,\!D_{1},\!D_{2}))&\mspace{-1.0mu}0\\ \gamma_{1}\lambda_{Z}(D_{2})f_{1}^{\prime}(N(\mu,\!D_{1},\!D_{2}))&\mspace{-15.0mu}\begin{matrix}\Bigl{(}\gamma_{1}f_{1}(N(\mu,\!D_{1},\!D_{2})){-}D_{1}\mspace{15.0mu}\\ \mspace{45.0mu}{-}Z(\mu,\!D_{1},\!D_{2})f_{2}^{\prime}(\lambda_{Z}(D_{2}))\Bigr{)}\end{matrix}&\mspace{-1.0mu}-D_{2}/\gamma_{2}\\ 0&\mspace{-15.0mu}\gamma_{2}Z(\mu,\!D_{1},\!D_{2})f_{2}^{\prime}(\lambda_{Z}(D_{2}))&\mspace{-1.0mu}0\end{bmatrix}.$$

$$x^3+a_1{x}^2+a_{2}x+a_3=0,$$

Full text

1 Introduction

We consider a mathematical model of two-species predator-prey interaction in the chemostat under nutrient limitation. With the exception of one nutrient, all nutrients that the prey species requires are supplied to the growth vessel from the feed vessel in ample supply. The predator species grows exclusively on the prey. With $N$ the concentration of the limiting nutrient, $P$ the concentration of prey (say, phytoplankton), and $Z$ the concentration of predator (say, zooplankton), we consider the following model:

[TABLE]

for initial conditions $N(0)=N_{0}>0$, $P(0)=P_{0}\geq 0$, and $Z(0)=Z_{0}\geq 0$.

The concentration of the growth-limiting nutrient in the feed vessel is denoted $\mu$, and will be the bifurcation parameter in our analysis. $D$ is the input rate from the feed vessel to the growth vessel as well as the washout rate from the growth vessel to the receptacle, so that constant volume is maintained. The parameters $D_{1}=D+\epsilon_{1}$ and $D_{2}=D+\epsilon_{2}$ are the removal rates of $P$ and $Z$, respectively, from the growth vessel, incorporating the washout rate $D$ and the intrinsic death rates $\epsilon_{i}$ of $P$ and $Z$. Our analysis does not necessarily require that $\epsilon_{1}$ and $\epsilon_{2}$ are positive; however, $D_{1}$ and $D_{2}$ should be positive. The yield coefficient $\gamma_{1}$ gives the amount of prey biomass produced per unit of nutrient consumed, while $\gamma_{2}$ gives the amount of predator biomass produced per unit of prey biomass consumed.

The function $f_{1}(N)$ represents the per capita consumption rate of nutrient by the prey populations as a function of the concentration of available nutrient; similarly, the function $f_{2}(P)$ represents the per capita consumption rate of the prey by the predator as a function of available prey. These functions are assumed to satisfy $f_{i}(0)=0$, $i=1$, 2, $f_{1}^{\prime}(N)>0$ for all $N\geq 0$, and $f_{2}^{\prime}(P)>0$ for all $P\geq 0$. We further assume that $f_{1}(N)$ and $f_{2}(P)$ are bounded. To avoid the case of washout due to an inadequate resource, we assume that $\lim_{N\to\infty}f_{1}(N)>D_{1}/\gamma_{1}$, and to avoid the case of an inadequate prey, we assume that $\lim_{P\to\infty}f_{2}(P)>D_{2}/\gamma_{2}$. Define $\lambda_{P}(D_{1})$ and $\lambda_{Z}(D_{2})$ to be the unique numbers satisfying

[TABLE]

The number $\lambda_{P}(D_{1})$ is the break-even concentration of nutrient, at which the growth and removal of phytoplankton balance; $\lambda_{Z}(D_{2})$ is similarly interpreted. The number $D$ plays a central role in our investigation, so we assume that $\lambda_{P}(D)$ and $\lambda_{Z}(D)$ are also defined. From the perspective of $D$, due to the boundedness assumptions on $f_{1}$ and $f_{2}$, $D_{1}$ and $D_{2}$ are perturbations of $D$.

Lemma 1.1**.**

From a functional point of view, $\lambda_{P}$ and $\lambda_{Z}$ are right inverses of $\gamma_{1}f_{1}$ and $\gamma_{2}f_{2}$, respectively. Accordingly, on their respective domains $\lambda_{P}$ and $\lambda_{Z}$ are as differentiable as $f_{1}$ and $f_{2}$. We note for later use

[TABLE]

Kuang and Li [7] studied this system with general functional responses and distinct values of $D$, $D_{1}$, and $D_{2}$. However, they fixed the input nutrient concentration, whereas we have this as a parameter. With the hypothesis that $D{=}D_{1}{=}D_{2}$, they provide numerical criteria for the stability of a coexistence equilibrium, and prove that a cycle exists when the equilibrium is unstable [7, Theorem 3.2]. When the hypothesis $D{=}D_{1}{=}D_{2}$ does not hold, they provide numerical evidence that stability of the coexistence equilibrium breaks down and a cycle appears. A similar model was studied in [9] with functional response $f_{1}(N)$ of Holling type I and $f_{2}(P)$ of Holling type II, demonstrating the existence of a Hopf bifurcation in response to varying nutrient concentration. These results inspire our work. Our goal is to prove analytically that the system undergoes a Hopf bifurcation without restricting the forms of the uptake functions or the values of the removal rates $D_{1}$ and $D_{2}$.

This paper is organized as follows. In section 2, conditions for existence and local stability of predator-free equilibria are obtained for general functional responses $f_{1}(N)$ and $f_{2}(P)$. In section 3, we study stability of a coexistence equilibrium which appears as the parameter $\mu$ increases. We quote a version of the Hopf bifurcation theorem, stating the result in a limited form most useful for us. Accordingly, application of this theorem requires us to develop specific information about the behavior of eigenvalues of the linearizations as the paramter $\mu$ increases. In section 4, we use a smooth change of variables to reach a situation where the Hopf bifurcation theorem in three dimensions applies, enabling us to conclude the existence of cycles. In section 5, our results are illustrated using simulations arising by choosing rate functions of Holling type II and Holling type III forms. In section 6, we explain in detail the approximations and estimates that support our work in the latter part of section 3.

2 Steady States and Their Stability

To begin, we establish that the solutions of (1) are nonnegative and bounded. These are minimum requirements for a reasonable model of the chemostat. We then develop conditions for the existence and stability of equilibria. We conclude the section by proving uniform persistence in the sense of [1] when $\mu$ is sufficiently large and the initial values $P_{0}$ and $Z_{0}$ are positive.

Lemma 2.1**.**

All solutions $N(t)$, $P(t)$, and $Z(t)$ of (1) are nonnegative and bounded.

Proof.

The plane $Z=0$ is invariant for system (1). Therefore, by existence and uniqueness, if $Z_{0}>0$ then $Z(t)>0$ for all $t\geq 0$. Similarly, since $f_{2}(0)=0$, the plane $P=0$ is invariant, so $P_{0}>0$ implies $P(t)>0$ for all $t\geq 0$.

Suppose $N_{0}>0$. If there exists a $t>0$ with $N(t)=0$, then there is a least such number, say $t_{0}$. Then $N^{\prime}(t_{0})=\mu D>0$ since $f_{1}(0)=0$. Consequently, there is $t<t_{0}$ such that $N(t)<0$, which is a contradiction to the choice of $t_{0}$.

For the boundedness of solutions, set $U(t)=N(t)+\gamma_{1}^{-1}P(t)+(\gamma_{1}\gamma_{2})^{-1}Z(t)$. From system (1), it follows that

[TABLE]

where $\widehat{D}=\min\{D,D_{1},D_{2}\}$. Then

[TABLE]

Since $N(t)$, $P(t)$, and $Z(t)$ are nonnegative, the boundedness of $U(t)$ implies the boundedness of $N(t)$, $P(t)$, and $Z(t)$. ∎

There are at most three biologically relevant equilibria of system (1) depending on the value of $\mu$. The equilibria and the conditions of their existence are summarized in the following theorem.

Theorem 2.2**.**

Let $\lambda_{P}(D_{1})$ and $\lambda_{Z}(D_{2})$ be as in (2). The equilibria of the system (1) satisfy the following conditions:

The washout equilibrium $E_{0}=(\mu,0,0)$ exists for all $\mu>0$. 2. 2.

The single-species equilibrium $E_{1}(\mu,D_{1})=(\lambda_{P}(D_{1}),P(\mu,D_{1}),0)$ exists for all $\mu>\lambda_{P}(D_{1})$, where

[TABLE] 3. 3.

The coexistence equilibrium

[TABLE]

exists for all $\displaystyle{\mu>\mu_{c_{1}}(D_{1},D_{2})}$, where

[TABLE]

and $N(\mu,D_{1},D_{2})$, $Z(\mu,D_{1},D_{2})$ satisfy the simultaneous equations

[TABLE]

Thus, for $\mu\leq\lambda_{P}(D_{1})$ only the equilibrium $E_{0}$ exists; for $\lambda_{P}(D_{1})<\mu\leq\mu_{c_{1}}(D_{1},D_{2})$, there are two equilibria $E_{0},E_{1}$; and, when $\mu_{c_{1}}(D_{1},D_{2})<\mu$, there are three equilibria $E_{0}$, $E_{1}$, and $E_{2}$.

Proof.

From the $Z$-equation, either $Z=0$ or $\gamma_{2}f_{2}(P)-D_{2}=0$ (so that $P=\lambda_{Z}(D_{2})$). If $Z=0$, the $P$-equation yields

[TABLE]

which implies either $P=0$ or $\gamma_{1}f_{1}(N)-D_{1}=0$ (so that $N=\lambda_{P}(D_{1})$).

If $Z=0$ and $P=0$, then the $N$-equation gives $0=(\mu-N)D$, so that $N=\mu$. Thus, the washout equilibrium is given by $E_{0}=(\mu,0,0)$ and exists for $\mu>0$. Note that as $\mu$ increases, $E_{0}$ moves along the $N$-axis in ${\mathbf{R}}_{\geq 0}^{3}$. This is the proof of the first part of the theorem.

If $Z=0$ and $P\neq 0$, then $N=\lambda_{P}(D_{1})$ in the $N$-equation gives

[TABLE]

Note that $P(\mu,D_{1})>0$ for all $\mu>\lambda_{P}(D_{1})$, and so the single-species equilibrium $E_{1}(\mu,D_{1})=(\lambda_{P}(D_{1}),P(\mu,D_{1}),0)$ exists in the positive cone for all $\mu>\lambda_{P}(D_{1})$. This is the proof of the second part of the theorem.

If $Z\neq 0$, then $P=\lambda_{Z}(D_{2})$ in the $N$- and $P$-equations, and equations (6) and (7) define $N$ and $Z$ as implicit functions of $\mu$, $D_{1}$, and $D_{2}$.

We now determine the critical value $\mu_{c_{1}}(D_{1},D_{2})$ of $\mu$ at which $E_{2}(\mu,D_{1},D_{2})$ first appears in the positive cone. Equation (7) implies

[TABLE]

Since $f_{1}$ is strictly increasing, $Z(\mu,D_{1},D_{2})=0$ if and only if $N(\mu,D_{1},D_{2})=\lambda_{P}(D_{1})$, and $Z(\mu,D_{1},D_{2})>0$ for $N(\mu,D_{1},D_{2})>\lambda_{P}(D_{1})$. Substituting $N(\mu,D_{1},D_{2})=\lambda_{P}(D_{1})$ in (6), we obtain the equation

[TABLE]

Thus, $Z(\mu,D_{1},D_{2})$ is positive when $\mu>\mu_{c_{1}}(D_{1},D_{2})$, and the coexistence equilibrium E_{2}(\mu,D_{1},D_{2})=\bigl{(}N(\mu,D_{1},D_{2}),\lambda_{Z}(D_{2}),Z(\mu,D_{1},D_{2})\bigr{)} exists in the positive cone for all $\mu>\mu_{c_{1}}(D_{1},D_{2})$. This proves part three of the theorem. ∎

In Theorems 2.3, 2.4, and 2.5 we investigate the stability of the equilibria of system (1) by finding the eigenvalues of the associated Jacobian matrices. The Jacobian matrix of the system (1) takes the form

[TABLE]

We first summarize the stability of $E_{0}$ in the following theorem. Here, the breakeven concentration of nutrient given in (2) plays a critical role.

Theorem 2.3**.**

The equilibrium point $E_{0}$ is locally asymptotically stable if $\mu<\lambda_{P}(D_{1})$ and unstable if $\mu>\lambda_{P}(D_{1})$. When $\mu>\lambda_{P}(D_{1})$, $E_{0}$ is globally asymptotically stable with respect to solutions initiating in $\{(N,P,Z)\in{\mathbf{R}}^{3}_{+}\mid P=0\}$. That is, the plane $P=0$ is $m^{+}(E_{0})$, the stable manifold of $E_{0}$.

Proof.

The Jacobian at $E_{0}=(\mu,0,0)$ is

[TABLE]

so that the eigenvalues of $J(E_{0})$ are $x_{1}=-D,x_{2}=\gamma_{1}f_{1}(\mu)-D_{1}$, and $\displaystyle{x_{3}=-D_{2}}$. The stability of $E_{0}$ now follows from (2) and the fact that $f_{1}$ is strictly increasing: $x_{2}<0$ when $\mu<\lambda_{P}(D_{1})$ and $x_{2}>0$ when $\mu>\lambda_{P}(D_{1})$. To see that $m^{+}(E_{0})=\{(N,P,Z)\in{\mathbf{R}}^{3}_{+}\mid P=0\}$ when $\mu>\lambda_{P}(D_{1})$, consider the Lyapunov function

[TABLE]

Clearly, $L(\mu,0)=0$ and $L(N,Z)>0$ if $(N,Z)\neq(\mu,0)$. The time derivative of $L(N,Z)$ at a point $(N,0,Z)$ on a trajectory of system (1) is

[TABLE]

for $(N,Z)\neq(\mu,0)$. Thus $E_{0}$ is globally asymptotically stable in the plane $P=0$. ∎

For $\mu=\lambda_{P}(D_{1})$, $P(\lambda_{P}(D_{1}),D_{1})=0$, so that $E_{0}$ and $E_{1}$ coalesce (see equation (4)). When $\mu>\lambda_{P}(D_{1})$, $E_{1}(\mu,D_{1})=(\lambda_{P}(D_{1}),P(\mu,D_{1}),0)$ enters the positive cone. We summarize the stability of $E_{1}(\mu,D_{1})$ in the following theorem. Note that the critical value of $\mu$ given in (5) now plays a central role.

Theorem 2.4**.**

The equilibrium point $E_{1}(\mu,D_{1})$ is locally stable if $\lambda_{P}(D_{1})<\mu<\mu_{c_{1}}(D_{1},D_{2})$ and unstable if $\mu>\mu_{c_{1}}(D_{1},D_{2})$. When $\mu>\mu_{c_{1}}(D_{1},D_{2})$, $E_{1}(\mu,D_{1})$ is globally asymptotically stable with respect to solutions initiating in $\{(N,P,Z)\in{\mathbf{R}}^{3}_{+}\mid Z=0\}$. That is, the plane $Z=0$ is m^{+}\bigl{(}E_{1}(\mu,D_{1})\bigr{)}, the stable manifold of $E_{1}(\mu,D_{1})$.

Proof.

The Jacobian matrix at $E_{1}(\mu,D_{1})$ is

[TABLE]

The determinant of the upper lefthand $2$-by-$2$ submatrix is positive and its trace is negative, so its eigenvalues have negative real parts. The third eigenvalue is x_{3}=\gamma_{2}f_{2}\bigl{(}P\bigl{(}\mu,D_{1})\bigr{)}-D_{2}. If $\mu<\mu_{c_{1}}(D_{1},D_{2})$, so that $P(\mu,D_{1})<\lambda_{Z}(D_{2})$, then $x_{3}<0$, and $E_{1}(\mu,D_{1})$ is locally stable. Similarly, if $\mu>\mu_{c_{1}}(D_{1},D_{2})$, so that $P(\mu,D_{1})>\lambda_{Z}(D_{2})$, then $x_{3}>0$ and $E_{1}(\mu,D_{1})$ is unstable with one dimension of instability.

To see that $m^{+}(E_{1})=\{(N,P,Z)\in{\mathbf{R}}^{3}_{+}\mid Z=0\}$ when $\mu>\mu_{c_{1}}$, consider the Lyapunov function introduced by Hsu in [4]

[TABLE]

Notice that L\bigl{(}\lambda_{P}(D_{1}),P(\mu,D_{1})\bigr{)}=0,

[TABLE]

precisely when (N,P)=\bigl{(}\lambda_{P}(D_{1}),P(\mu,D_{1})\bigr{)}. Moreover,

[TABLE]

Therefore, \bigl{(}\lambda_{P}(D_{1}),P(\mu,D_{1})\bigr{)} is the only critical point of $L(P,N)$ and it is a local minimum, so that $L(N,P)>0$ for all (N,P)\neq\bigl{(}\lambda_{P}(D_{1}),P(\mu,D_{1})\bigr{)}.

Now we compute the time derivative of $L(N,P)$ at a point $(N,P,0)$ along a trajectory of system (1). Noting from (2) and (4) that

[TABLE]

we have

[TABLE]

If $N<\lambda_{P}(D_{1})$, then f_{1}(N)-f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}<0, so that f_{1}(N)/f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}<1. Also, $\mu{-}N>\mu{-}\lambda_{P}(D_{1})$, so that

[TABLE]

Thus, $L^{\prime}(N,P)<0$ when $N<\lambda_{P}(D_{1})$.

If $N>\lambda_{P}(D_{1})$, then f_{1}(N)-f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}>0, so that f_{1}(N)/f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}>1. When $\mu>N>\lambda_{P}(D_{1})$, we have $\mu{-}N<\mu{-}\lambda_{P}(D_{1})$, so that

[TABLE]

while for $N\geq\mu$, then $\displaystyle{\frac{\mu-N}{\mu-\lambda_{P}(D_{1})}<0}$. In either case, \displaystyle{\frac{\mu-N}{\mu-\lambda_{P}(D_{1})}-\frac{f_{1}(N)}{f_{1}\bigl{(}\lambda_{P}(D_{1})\bigr{)}}<0}, and so $L^{\prime}(N,P)<0$ when $N>\lambda_{P}(D_{1})$.

Finally, $L^{\prime}(N,P)=0$ if and only if $N=\lambda_{P}(D_{1})$. By LaSalle’s extension theorem [6], any trajectory of system (1) in the plane $Z=0$ for which $P_{0}>0$ approaches the largest invariant set in the line $N=\lambda_{P}(D_{1})$, and this is simply $\{E_{1}(\mu,D_{1})\}$. Therefore, $E_{1}(\mu,D_{1})$ is globally asymptotically stable in the plane $Z=0$. ∎

When $\mu=\mu_{c_{1}}(D_{1},D_{2})$, $E_{1}(\mu,D_{1})$ and $E_{2}(\mu,D_{1},D_{2})$ coalesce. With $\mu=\mu_{c_{1}}(D_{1},D_{2})$, we have P\bigl{(}\mu_{c_{1}}(D_{1},D_{2}),D_{1}\bigr{)}=\lambda_{Z}(D_{2}) (see equations (5) and (4)). Also $N(\mu,D_{1},D_{2})=\lambda_{P}(D_{1})$, so that $Z(\mu,D_{1},D_{2})=0$ (see the discussion around (9)). Thus,

[TABLE]

Said another way, as $\mu$ increases through $\mu_{c_{1}}(D_{1},D_{2})$, $E_{2}(\mu,D_{1},D_{2})$ enters the positive cone by passing through $E_{1}(\mu,D_{1})$.

With f_{2}\bigl{(}\lambda_{Z}(D_{2})\bigr{)}=D_{2}/\gamma_{2}, the Jacobian matrix at $E_{2}(\mu,D_{1},D_{2})$ takes the form

[TABLE]

The eigenvalues of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} satisfy

[TABLE]

where

[TABLE]

Theorem 2.5**.**

The coexistence equilibrium point

[TABLE]

is asymptotically stable if and only if $a_{1}>0$ and $a_{1}a_{2}>a_{3}$.

Proof.

Since $a_{3}$ is positive, this follows from the Routh-Hurwitz criterion. ∎

We conclude this section with a significant strengthening of Lemma 2.1. We use the concept of uniform persistence introduced in [1].

Theorem 2.6**.**

Assume that $\mu>\mu_{c_{1}}(D_{1},D_{2})$. Then system (1) is uniformly persistent with respect to all solutions satisfying $P_{0}>0$ and $Z_{0}>0$.

Proof.

Recall from Lemma 2.1 that all solutions of system (1) for which $P_{0}>0$ and $Z_{0}>0$ are positive and bounded.

We first show that $\liminf_{t\rightarrow\infty}{N(t)>0}$. If $\liminf_{t\rightarrow\infty}N(t)=0$ and $\limsup_{t\rightarrow\infty}N(t)=0$, then $\lim_{t\rightarrow\infty}N(t)=0$. But this is impossible, for then it follows from the $N$-equation that $N^{\prime}(t)\rightarrow\mu D>0$ as $t\rightarrow\infty$ and this, in turn, contradicts the fact that $N(t)$ is bounded.

Now, suppose $\liminf_{t\rightarrow\infty}{N(t)=0}$ while $\limsup_{t\rightarrow\infty}{N(t)>0}$. Then there exists a sequence $\{\tau_{n}\}_{n=1}^{\infty}$ of local minima of $N(t)$ satisfying $\tau_{n}\rightarrow\infty$ as $n\rightarrow\infty$. Thus,

$N^{\prime}(\tau_{n})=0$, since $\tau_{n}$ is a local minimum, and 2. 2.

$N(\tau_{n})\rightarrow 0$ as $n\rightarrow\infty$, since $\liminf_{t\rightarrow\infty}{N(t)=0}$.

From the $N$-equation we have

[TABLE]

so that

[TABLE]

Rearranging and using the facts that $N(\tau_{n})\rightarrow 0$ as $n\rightarrow\infty$, $f_{1}$ is continuous and $f_{1}(0)=0$, we get

[TABLE]

a contradiction. Hence, $\liminf_{t\rightarrow\infty}{N(t)>0}$.

Choose X(0)=\bigl{(}N_{0},P_{0},Z_{0}\bigr{)}\in{\mathbf{R}}^{3}_{+}. Then \omega\bigl{(}X(0)\bigr{)} is a nonempty, compact invariant set with respect to system (1). We claim $E_{0}=(\mu,0,0)$ and $E_{1}(\mu,D_{1})=(\lambda_{P}(D_{1}),P(\mu,D_{1}),0)$ are not in $\omega(X(0))$. Suppose E_{0}=(\mu,0,0)\in\omega\bigl{(}X(0)\bigr{)}. Since $\mu>\mu_{c_{1}}(D_{1},D_{2})$, $E_{0}$ is an unstable hyperbolic equilibrium point. By theorem 2.3, $E_{0}$ is globally asymptotically stable with respect to solutions initiating in the plane $P=0$. Since $X(0)\not\in m^{+}(E_{0})$, \{E_{0}\}\neq\omega\bigl{(}X(0)\bigr{)}. By the Butler-McGehee lemma (see lemma A1, [2]), there exists Q\in\bigl{(}m^{+}(E_{0})\setminus\{E_{0}\}\bigr{)}\cap\omega\bigl{(}X(0)\bigr{)}, so that $\operatorname{cl}\mathcal{O}(Q)\subset\omega(X(0))$. For such an initial condition, the governing system is \displaystyle{N^{\prime}(t)=D\bigl{(}\mu-N(t)\bigr{)}}, and $\displaystyle{Z^{\prime}(t)=-D_{2}Z(t)}$. But then $\mathcal{O}(Q)$ becomes unbounded as $t\rightarrow-\infty$. This is a contradiction to the compactness of \omega\bigl{(}X(0)\bigr{)}, and so E_{0}\not\in\omega\bigl{(}X(0)\bigr{)}.

Suppose $E_{1}(\mu,D_{1})=(\lambda_{P}(D_{1}),P(\mu,D_{1}),0)$ is in $\omega(X(0))$. Since $\mu>\mu_{c_{1}}(D_{1},D_{2})$, $E_{1}(\mu,D_{1})$ is an unstable hyperbolic equilibrium point. By theorem 2.4 $E_{1}(\mu,D_{1})$ is globally asymptotically stable with respect to solutions initiating in the plane $Z=0$. Since X(0)\not\in m^{+}\bigl{(}E_{1}(\mu,D_{1})\bigr{)}, \{E_{1}(\mu,D_{1})\}\neq\omega\bigl{(}X(0)\bigr{)}. By the Butler-McGehee lemma, there exists \widehat{Q}\in\bigl{(}m^{+}\bigl{(}E_{1}(\mu,D_{1})\bigr{)}\setminus\{E_{1}(\mu,D_{1})\}\bigr{)}\cap\omega\bigl{(}X(0)\bigr{)}, so that $\operatorname{cl}\mathcal{O}(\widehat{Q})\subset\omega(X(0))$. If $\widehat{Q}\in m^{-}(E_{0})$, then $E_{0}\in\operatorname{cl}{\mathcal{O}(\widehat{Q})}\subset\omega(X(0))$, a contradiction. Thus $\widehat{Q}\not\in m^{-}(E_{0})$, and this implies $\mathcal{O}(\widehat{Q})$ is unbounded as $t\rightarrow-\infty$, contradicting the compactness of $\omega(X(0))$. Therefore, $E_{1}(\mu,D_{1})\not\in\omega(X(0))$.

Suppose the system (1) is not persistent. Then there exists $\widetilde{Q}\in\omega(X(0))$ such that $\widetilde{Q}\in m^{+}(E_{0})$ or \widetilde{Q}\in m^{+}\bigl{(}E_{1}(\mu,D_{1})\bigr{)}. Then $\operatorname{cl}{\mathcal{O}(\widetilde{Q})}\subset\omega(X(0))$, which implies either $E_{0}\in\omega(X(0))$ or $E_{1}(\mu,D_{1})\in\omega(X(0))$, neither of which can be true. Thus, $\liminf_{t\rightarrow\infty}P(t)>0$ and $\liminf_{t\rightarrow\infty}Z(t)>0$, and it follows from the main result of [1] that system (1) is uniformly persistent. ∎

3 Stability of the coexistence equilibrium

In this section we study the coexistence equilibrium point

[TABLE]

first when $D{=}D_{1}{=}D_{2}$ and $\mu$ varies and then after relaxing the assumption $D{=}D_{1}{=}D_{2}$. In particular, we study the eigenvalues of the Jacobians at the coexistence equilibria as $\mu$ increases. To prepare for the study of the evolution of the equilibrium $E_{2}(\mu,D_{1},D_{2})$, we observe the following consequence of the implicit function theorem [5, p.122].

Lemma 3.1**.**

For $D_{1}>0$, $D_{2}>0$, and $\mu>\mu_{c_{1}}(D_{1},D_{2})$ there is a local parameterization of the locus of coexistence equilibria

[TABLE]

defined on an interval containing $\mu$ and a disc around $(D_{1},D_{2})$ that is smooth to the smaller of the degree of smoothness of $f_{1}(N)$ and the degree of smoothness of $f_{2}(P)$.

Proof.

Set $P=\lambda_{Z}(D_{2})$ in the $N$- and $P$- equations of system (1). With $f_{2}(\lambda_{Z}(D_{2}))=D_{2}/\gamma_{2}$ from (2), let

[TABLE]

and

[TABLE]

and define $G\colon{\mathbf{R}}^{5}\rightarrow{\mathbf{R}}^{2}$ by G(N,Z,\mu,D_{1},D_{2})=\bigl{(}G_{1}(N,Z,\mu,D_{1},D_{2}),G_{2}(N,Z,\mu,D_{1},D_{2})\bigr{)}. Then we want to parametrize the set $G^{-1}(0,0)$.

The derivative of $G$ is represented by the matrix

[TABLE]

Observe that, for fixed $\mu_{0}>\mu_{c_{1}}(D_{1},D_{2})$, the first two columns of $DG$ at the point

[TABLE]

evaluate to

[TABLE]

Since $f_{1}^{\prime}(N)>0$, the first two columns are linearly independent and so the implicit function theorem applies. There exists a ball $B_{1}$ around $(\mu_{0},D_{1},D_{2})$ and a ball $B_{2}$ around \bigl{(}N(\mu_{0},D_{1},D_{2}),Z(\mu_{0},D_{1},D_{2})\bigr{)} such that for each $(\mu,D_{1},D_{2})$ in $B_{1}$ there is a unique point \bigl{(}N(\mu,D_{1},D_{2}),Z(\mu,D_{1},D_{2})\bigr{)} in $B_{2}$ such that G\bigl{(}N(\mu,D_{1},D_{2}),Z(\mu,D_{1},D_{2}),\mu\bigr{)}=0. Moreover, the functions $N(\mu,D_{1},D_{2})$ and $Z(\mu,D_{1},D_{2})$ have the same degree of differentiability as does $G$, which is the minimum of the degrees of differentiability of $f_{1}(N)$ and $f_{2}(P)$. To explain how the differentiability of $f_{2}$ enters, note that the computation of $\partial N/\partial D_{2}$ and $\partial Z/\partial D_{2}$ via $\partial G/\partial D_{2}$ involves \lambda_{Z}^{\prime}(D_{2})=\bigl{(}\gamma_{2}\cdot f_{2}^{\prime}(\lambda_{P}(D_{2}))\bigr{)}^{-1} by (3). This all gives us a smooth parameterization of the equilibrium locus, as desired. ∎

Remark*.*

From the point of view of calculus, each of the independent variables $\mu$, $D_{1}$, $D_{2}$ is on an equal footing with the others. However, viewed through the lens of the structure of system (1), we can describe the variable $\mu$ as a control parameter and the variables $D_{1}$ and $D_{2}$ as experimental parameters fixed at some earlier time. This distinction informs our analysis where we first consider the model when $D{=}D_{1}{=}D_{2}$ as $\mu$ varies, and subsequently vary $D_{1}$ and $D_{2}$ in a neighborhood of $D$.

Theorem 3.2**.**

Fix $D_{1}$ and $D_{2}$ and suppose $\mu>\mu_{c_{1}}(D_{1},D_{2})$, so that the coexistence equilibrium point E_{2}(\mu,D_{1},D_{2})=\bigl{(}N(\mu,D_{1},D_{2}),\lambda_{Z}(D_{2}),Z(\mu,D_{1},D_{2})\bigr{)} exists. Then

the $\mu$-derivatives $N^{\prime}(\mu,D_{1},D_{2})$ and $Z^{\prime}(\mu,D_{1},D_{2})$ are positive. 2. 2.

$\lim_{\mu\rightarrow\infty}N(\mu,D_{1},D_{2})=\infty$* and $Z(\mu,D_{1},D_{2})$ is bounded.*

Proof.

Since $D_{1}$ and $D_{2}$ are fixed throughout this proof, we drop these symbols and write simply $N(\mu)$ and $Z(\mu)$. To show $N^{\prime}(\mu)>0$, recall equation (6) from theorem 2.2:

[TABLE]

Differentiating (6) with respect to $\mu$ and rearranging, we get

[TABLE]

Since $f_{1}^{\prime}(N)$ is positive, it follows that $N^{\prime}(\mu)$ is positive.

To show $Z^{\prime}(\mu)>0$, set f_{2}\bigl{(}\lambda_{Z}(D_{2})\bigr{)}=D_{2}/\gamma_{2} in equation (7), obtaining

[TABLE]

Differentiating with respect to $\mu$ and rearranging, we get

[TABLE]

Since $N^{\prime}(\mu)$ and $f_{1}^{\prime}(N)$ are positive, it follows that $Z^{\prime}(\mu)$ is positive. This proves part one.

To prove part two, note that equation (6) implies

[TABLE]

Since $f_{1}(N)$ is bounded, $\lim_{\mu\rightarrow\infty}N(\mu)=\infty$. From (7) we have

[TABLE]

Since $f_{1}(N)$ is bounded, $Z(\mu)$ is bounded. ∎

Before we turn to a study of the stability properties of the coexistence equilibrium, we include a partial paraphrase of the $C^{L}$ Hopf bifurcation theorem as stated in [3, p.16]. Since our goal is to make an application of this result to the coexistence equilibrium $E_{2}$, and because the verification of the hypotheses is lengthy, we refer to our paraphrase to keep track of progress.

Theorem 3.3**.**

Consider a system $dX/dt=F(X,\mu)$ with $X\in{\mathbf{R}}^{n}$ and $\mu$ a real parameter. If,

for $\mu$ in an open interval containing $\mu_{c}$ (characterized in 3 below), $F(0,\mu)=0$ and $0\in{\mathbf{R}}^{n}$ is an isolated equilibrium point of $dX/dt{=}F(X,\mu)$; 2. 2.

all partial derivatives of the components $F^{\ell}$ of the vector $F$ of orders $\leq L{+}2$, ($L\geq 2$) exist and are continuous in $X$ and $\mu$ in a neighborhood of $(0,\mu_{c})$ in ${\mathbf{R}}^{n}{\times}{\mathbf{R}}$; 3. 3.

the Jacobian $J(0,\mu)=D_{X}F(0,\mu)$ has a pair of complex eigenvalues $\alpha(\mu){\pm}i\,\omega(\mu)$, where $\alpha(\mu_{c})=0$ and $\alpha^{\prime}(\mu_{c})\neq 0$; 4. 4.

the remaining $n{-}2$ eigenvalues of $J(0,\mu_{c})$ have strictly negative real parts,

then the system $dX/dt=F(X,\mu)$ has a family of periodic solutions.

Remark*.*

For the purposes of the proof in [3] the authors assume the critical value of the bifurcation parameter is $\mu_{c}=0$, which is a trivial alteration. We are only interested in the existence of cycles, so we do not quote several additional conclusions offered in [3]. In our situation the equilibrium $E_{2}$ depends on the parameter $\mu$, an issue which we circumvent in section 4 using the inverse function theorem. Hypothesis 2 is fulfilled by imposing more differentiability conditions on the functions $f_{1}(N)$ and $f_{2}(P)$ at an appropriate point in the exposition. Verification of hypotheses 3 and 4 in the statement of theorem 3.3 is the most involved part of our process and occupies the rest of this section.

To begin the stability analysis, we conjugate the Jacobian $J(E_{2})$ in (11) by

[TABLE]

Using f_{2}\bigl{(}\lambda_{Z}(D_{2})\bigr{)}=D_{2}/\gamma_{2} and writing $D_{1}=D{+}\epsilon_{1}$, $D_{2}=D{+}\epsilon_{2}$ yields the matrix

[TABLE]

where

[TABLE]

First we make the assumption that $D{=}D_{1}{=}D_{2}$, so that $\epsilon_{1}{=}\epsilon_{2}{=}0$; that is, we assume the death rates of $P$ and $Z$ are negligible with respect to washout rate $D$. The results of [7] suggest this is a useful initial assumption. Since we regard $D$ as fixed for the discussion, we will abbreviate $N(\mu,D,D)$ by $N(\mu)$ and $Z(\mu,D,D)$ by $Z(\mu)$. Similarly, we will abbreviate $E_{2}(\mu,D,D)$ by $E_{2}(\mu)$.

We can explicitly compute the eigenvalues of J\bigl{(}E_{2}(\mu)\bigr{)}, since conjugation does not change them. The characteristic polynomial of J\bigl{(}E_{2}(\mu)\bigr{)} is

[TABLE]

where

[TABLE]

The next result amplifies theorem 2.5 in the case $D{=}D_{1}{=}D_{2}$.

Theorem 3.4**.**

Assume $D{=}D_{1}{=}D_{2}$ and that $\mu>\mu_{c_{1}}(D,D)$, so the coexistence equilibrium $E_{2}(\mu)=(N(\mu),\lambda_{Z}(D),Z(\mu))$ exists. Then

$E_{2}(\mu)$* is locally asymptotically stable if*

[TABLE] 2. 2.

$E_{2}(\mu)$* is unstable if*

[TABLE]

Proof.

With f_{2}\bigl{(}\lambda_{Z}(D)\bigr{)}=D/\gamma_{2} in (7) we have \gamma_{1}f_{1}\bigl{(}N(\mu)\bigr{)}-D=\bigl{(}DZ(\mu)\bigr{)}/\bigl{(}\gamma_{2}\lambda_{Z}(D)\bigr{)}. Then

[TABLE]

The result now follows from the Routh-Hurwitz criterion, since $-BC>0$ is easily verified from formulas (18) and (19). ∎

From the factorization of the characteristic polynomial of J\bigl{(}E_{2}(\mu)\bigr{)} given in (16) it is immediate that the Jacobian at $E_{2}(\mu)$ has one negative eigenvalue. Thus, hypothesis 4 of theorem 3.3 for $E_{2}(\mu)$ is satisfied in the case $D{=}D_{1}{=}D_{2}$. Now we verify hypothesis 3 for this situation.

Theorem 3.5**.**

Assume $D{=}D_{1}{=}D_{2}$ and let $\mu>\mu_{c_{1}}(D,D)$, so that the coexistence equilibrium E_{2}(\mu)=\bigl{(}N(\mu),\lambda_{Z},Z(\mu)\bigr{)} exists.

If $f^{\prime}_{1}(N)$ is continuous and D/\bigl{(}\gamma_{2}\lambda_{Z}(D)\bigr{)}>f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)}, then there exists a value $\mu_{c_{2}}>\mu_{c_{1}}(D,D)$ for which $A(\mu_{c_{2}})=0$. Consequently, when $\mu=\mu_{c_{2}}$, the Jacobian has a conjugate pair of imaginary eigenvalues. 2. 2.

If, in addition, $f_{1}$ is twice differentiable with respect to $N$ and f_{1}^{(2)}\bigl{(}N(\mu_{c_{2}})\bigr{)}<0, then $A^{\prime}(\mu_{c_{2}})>0$. Combining this with part 1, we have that hypothesis 3 of theorem 3.3 is satisfied for $E_{2}(\mu)$. 3. 3.

If $f_{1}^{(2)}(N)<0$ for all $N$, then $\mu_{c_{2}}$ is unique.

Remark*.*

If we assume $f^{(2)}_{2}(P)<0$, i.e., that $f_{2}$ is is concave down, then the slope of the secant that passes through the points $(0,0)$ and \bigl{(}\lambda_{Z}(D),(D/\gamma_{2})\bigr{)} is greater than the slope of the tangent line to the graph of $f_{2}$ at $\lambda_{Z}(D)$; that is, D/\bigl{(}\gamma_{2}\lambda_{Z}(D)\bigr{)}>f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)}. This will be the case, for example, when $f_{2}(P)$ is Holling Type II. However, the one-point condition D/\bigl{(}\gamma_{2}\lambda_{Z}(D)\bigr{)}-f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)}>0 may hold even if $f_{2}$ is not concave down. We will give such an example in Section 5.

Proof.

Consider the expression

[TABLE]

given in (20). Since $f_{1}^{\prime}(N)$ is a continuous function, $A$ is a continuous function of $\mu$ by Lemma 3.1. To prove that $A(\mu)$ has a zero value for some $\mu>\mu_{c_{1}}(D,D)$, it is enough to prove that $A(\mu)$ passes from negative to positive. For $\mu=\mu_{c_{1}}(D,D)$, Z\bigl{(}\mu_{c_{1}}(D,D)\bigr{)}=0, so

[TABLE]

To find a value of $\mu>\mu_{c_{1}}(D,D)$ at which $A(\mu)$ is positive, we use Theorem 3.2. We have that $Z(\mu)$ is increasing and bounded for $\mu>\mu_{c_{1}}(D,D)$. Let ${Z_{\infty}=\sup_{\mu\geq\mu_{c_{1}}(D,D)}Z(\mu)}$. Then there exists an $M_{1}$ such that for $\mu>M_{1}$, $Z(\mu)>Z_{\infty}/2$. Since $f_{1}$ is bounded and increasing, $\lim_{N\to+\infty}f^{\prime}_{1}(N)=0$. Then, for any $\epsilon>0$, there exists an $N_{\epsilon}>0$ such that $0<f_{1}^{\prime}(N)<\epsilon$ for all $N>N_{\epsilon}$. With \epsilon=\bigl{(}Z_{\infty}/2\lambda_{Z}(D)\bigr{)}\cdot\Bigl{(}D/\bigl{(}\gamma_{2}\lambda_{Z}(D)\bigr{)}-f^{\prime}_{2}\bigl{(}\lambda_{Z}(D)\bigr{)}\Bigr{)}, this implies that there exists an $N^{*}$ such that, if $N>N^{*}$, then

[TABLE]

In addition, $N(\mu)$ is increasing without bound by theorem 3.2, so there is an $M_{2}>\mu_{c_{1}}(D,D)$ such that, if $\mu>M_{2}$, then $N(\mu)>N^{*}$. Choose $\mu^{*}>\max\{M_{1},M_{2}\}$. Then

[TABLE]

Since A\bigl{(}\mu_{c_{1}}(D,D)\bigr{)}<0 and $A(\mu^{*})>0$, there is a number $\mu_{c_{2}}>\mu_{c_{1}}(D,D)$ such that $A(\mu_{c_{2}})=0$. Note that when $\mu=\mu_{c_{2}}$ the discriminant of the quadratic factor of the characteristic polynomial in (16) is

[TABLE]

so its roots are indeed purely imaginary. This proves part one.

For part two, by continuity of the discriminant as a function of $\mu$, there is a neighborhood of $\mu_{c_{2}}$ on which the discriminant is negative. Continuing, differentiate $A(\mu)$ with respect to $\mu$ to obtain

[TABLE]

By theorem 3.2, $N^{\prime}(\mu)$ and $Z^{\prime}(\mu)$ are positive. By the hypotheses of the present theorem, D/\bigl{(}\gamma_{2}\lambda_{Z}(D)\bigr{)}-f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)}>0 and $f_{1}^{(2)}(N(\mu_{c_{2}}))<0$. Thus, $A^{\prime}(\mu_{c_{2}})>0$. Combining parts one and two means that hypothesis 3 of Theorem 3.3 holds at $\mu_{c_{2}}$ for the case $D{=}D_{1}{=}D_{2}$.

For part three, if $f_{1}^{(2)}(N)<0$ for all $N$, then $A^{\prime}(\mu)>0$ for $\mu>\mu_{c_{1}}(D,D)$. ∎

Let us now discuss weakening the condition $D{=}D_{1}{=}D_{2}$. Intuitively, for $(D_{1},D_{2})$ sufficiently close to $(D,D)$ the eigenvalues of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} should exhibit behavior similar to those of J\bigl{(}E_{2}(\mu,D,D)\bigr{)}. In particular, the equilibrium $E_{2}(\mu,D_{1},D_{2})$ should exhibit a similar loss of stability. To make these considerations precise, we first require lemma 3.6.

Lemma 3.6**.**

Let $P_{1}(x)^{-}=\{(\alpha-x)\mid\alpha<0\}$ be the space of polynomials of degree 1 in $x$, with leading coefficient $-1$ and negative constant term, let $P_{2}(x)^{-}=\{\beta-\gamma x+x^{2}\mid\gamma^{2}-4\beta<0\}$ be the space of monic quadratic polynomials in $x$ with real coefficients and having a complex conjugate pair of roots, and let $P_{3}(x)^{-}=\{p_{0}+p_{1}x+p_{2}x^{2}-x^{3}\}$ be the space of cubic polynomials in $x$ with leading coefficient $-1$ and real coefficients.

Then the multiplication map $M\colon P_{1}^{-}\times P_{2}^{-}\rightarrow P_{3}^{-}$ is locally a diffeomorphism.

Proof.

Identify $P_{3}^{-}$ with Euclidean space using $p_{0}$, $p_{1}$, and $p_{2}$ as coordinates, identify $P_{1}^{-}$ with an open subset of ${\mathbf{R}}$ using $\alpha$ as the coordinate, and identify $P_{2}^{-}$ with an open subset of the plane ${\mathbf{R}}^{2}$ using $\beta$ and $\gamma$ as coordinates. Then the map $M$ has the expression

[TABLE]

The derivative, or Jacobian, of $M$ at $(\alpha,\beta,\gamma)$ is

[TABLE]

which fails to be invertible if and only if

[TABLE]

Should this occur, then \alpha=\bigl{(}\gamma\pm\sqrt{\gamma^{2}-4\beta}\bigr{)}/2. But we assume $\alpha<0$ is real and $\gamma^{2}-4\beta<0$, so $\det(DM)=0$ is impossible. The map $M$ is smooth, so, by the inverse function theorem [5, p.125], it is locally a diffeomorphism. ∎

To explain how lemma 3.6 comes into play, consider the map $\chi\colon M_{3,3}({\mathbf{R}})\rightarrow P_{3}^{-}$ which takes as input a real-valued $3$ by $3$ matrix $R$ and produces its characteristic polynomial $\chi(R)=\det(R-xI)$. The map $\chi$ has a coordinate expression by taking the coefficients in degrees [math], $1$, and $2$. These coefficients are polynomials in the matrix entries, so $\chi$ is smooth. Now look at

[TABLE]

Suppose we are in the situation of theorem 3.5, where it is easy to factor the characteristic polynomial of J\bigl{(}E_{2}(\mu,D,D)\bigr{)}, and we have seen the Jacobian has a negative real eigenvalue and a complex conjugate pair of eigenvalues as $\mu$ varies near a potential bifurcation value $\mu_{c_{2}}$. The factorization explicitly inverts the polynomial multiplication $M$ at particular points, and lemma 3.6 enables us to extend the factorization, in principle, to the characteristic polynomial of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)}.

In particular, we can understand how the characteristic polynomial of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} behaves as $\mu$ varies when $(D_{1},D_{2})$ is close to $(D,D)$ (in the Euclidean norm, for definiteness). We remind the reader that we think of $\mu$ as a control parameter, adjustable by the experimenter, and $D_{1}$ and $D_{2}$ as experimental parameters, set at the beginning of an experiment. To bring out this distinction, we will write the components of the formal factorization of the characteristic polynomial of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} as $\alpha(D_{1},D_{2})(\mu)$, $\beta(D_{1},D_{2})(\mu)$, and $\gamma(D_{1},D_{2})(\mu)$.

Lemma 3.7**.**

Assume $f_{1}$ is three times continuously differentiable and $f_{2}$ is two times continuously differentiable. Then there exists a $\mu$-interval $[\mu{-}\hat{\delta},\mu{+}\hat{\delta}]$ on which the $\mu$-derivative $\gamma^{\prime}(D_{1},D_{2})(\mu)$ is uniformly approximated by $\gamma^{\prime}(D,D)(\mu)=A^{\prime}(\mu)$. In fact, there exists a constant $C$ such that

[TABLE]

for any $\mu\in[\mu_{c_{2}}{-}\hat{\delta},\mu_{c_{2}}{+}\hat{\delta}]$

The details of the proof of lemma 3.7 are relegated to section 6 so as not to disturb the flow of the exposition.

Theorem 3.8**.**

Let $\mu>\mu_{c_{1}}(D_{1},D_{2})$, so that the equilibrium point $E_{2}(\mu,D_{1},D_{2})$ exists in the interior of the positive octant. Assume $f_{1}$ is three times continuously differentiable and $f_{2}$ is two times continuously differentiable and that D/\bigl{(}\gamma_{2}\lambda_{Z}(D)\bigr{)}>f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)} so that the hypotheses of theorem 3.5 part 1 are satisfied for $E_{2}(\mu,D,D)$, and let $\mu_{c_{2}}$ be as in theorem 3.5.

For $(D_{1},D_{2})$ sufficiently close to $(D,D)$,

hypothesis 4 of theorem 3.3 holds for $E_{2}(\mu,D_{1},D_{2})$; 2. 2.

and if, in addition, f_{1}^{(2)}\bigl{(}N(\mu_{c_{2}},D,D)\bigr{)}<0 so that the hypotheses of theorem 3.5 part 2 are satisfied, then hypothesis 3 of theorem 3.3 holds for $E_{2}(\mu,D_{1},D_{2})$.

Remark*.*

Theorem 3.5 gives a condition, namely, $f_{1}^{(2)}<0$, on system (1) under which there is a unique number $\mu_{c_{2}}$ at which J\bigl{(}E_{2}(\mu,D,D)\bigr{)} has a purely imaginary pair of eigenvalues meeting the transverality condition. However, it is not a priori clear that, in general, there is precisely one number at which these properties of the eigenvalues hold. Therefore, for the theorem and its proof, choose one such number $\mu_{c_{2}}$ and fix it throughout the discussion.

Proof of theorem 3.8..

We have proved in theorem 3.5 that there is an interval of parameters $\mu$ in which the characteristic polynomials

[TABLE]

of J\bigl{(}E_{2}(\mu,D,D)\bigr{)} have a complex conjugate pair of roots in addition to the eigenvalue $-D{<}0$, so they are in the image of the multiplication map $M$. For $(D_{1},D_{2})$ close to $(D,D)$ the entries in J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} are close to the entries in J\bigl{(}E_{2}(\mu,D,D)\bigr{)}, so the characteristic polynomial of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} is close to the characteristic polynomial of $J(E_{2})(\mu,D,D)$. To see this explicitly, refer to the formulas (12), (13), and (14); to account for the normalization of the characteristic polynomial to leading coefficient $-1$ multiply each expression by $-1$. Therefore, in view of lemma 3.6, the characteristic polynomial of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} has a decomposition of the same form as that of the characteristic polynomial of J\bigl{(}E_{2}(\mu,D,D)\bigr{)}. Written formally, the decomposition is

[TABLE]

The map $M\colon P_{1}^{-}(x){\times}P_{2}^{-}(x)\rightarrow P_{3}^{-}(x)$ is infinitely differentiable, so the local inverse $M^{-1}$ is also. In particular, for a fixed value $\mu$, the coefficients of the decomposition are smooth functions of $(D_{1},D_{2})$ defined in a neighborhood of $(D,D)$.

The first consequence is that, by definition of $P_{1}^{-}$, the characteristic polynomial of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} has a linear factor $\alpha(D_{1},D_{2})(\mu)-x$ with $\alpha(D_{1},D_{2})(\mu)<0$. This shows that hypothesis 4 of theorem 3.3 can be satisfied.

Now we start the verification of hypothesis 3 of theorem 3.3 for $E_{2}(\mu,D_{1},D_{2})$. The discriminant of the quadratic factor of the characteristic polynomial of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} is

[TABLE]

At $\mu_{c_{2}}$ and for $(D_{1},D_{2})$ sufficiently close to $(D,D)$, this is close to the expression

[TABLE]

for the discriminant of the quadratic factor of the characteristic polynomial of J\bigl{(}E_{2}(\mu,D,D)\bigr{)} given in (21), which is negative. Therefore, at $\mu_{c_{2}}$, the discriminant \bigl{(}\gamma(D_{1},D_{2})(\mu)\bigr{)}^{2}-4\,\beta(D_{1},D_{2})(\mu) is also negative. By continuity of the discriminant as a function of $\mu$, there is an interval $[\mu_{c_{2}}-\delta_{0},\mu_{c_{2}}+\delta_{0}]$ on which it is negative. Therefore, the characteristic polynomial of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} has a complex conjugate pair of roots on this interval.

Under the assumptions that $f_{1}^{(2)}$ is continuous and f_{1}^{(2)}\bigl{(}N(\mu_{c_{2}},D,D)\bigr{)}<0, $A^{\prime}(\mu)$ as calculated in (22) is continuous and $A^{\prime}(\mu_{c_{2}})>0$. So there is a $\delta_{1}>0$ such that $A^{\prime}(\mu)>A^{\prime}(\mu_{c_{2}})/2$ on the interval $[\mu_{c_{2}}-\delta_{1},\mu_{c_{2}}+\delta_{1}]$. Choose $\delta=\min\{\delta_{0},\delta_{1}\}$. Put $\eta_{1}=-A(\mu_{c_{2}}{-}\delta)/2>0$. By continuity of $\gamma$ as a function of $(D_{1},D_{2})$, there is a $\rho_{1}>0$ such that {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}<\rho_{1} implies

[TABLE]

Similarly, put $\eta_{2}=A(\mu_{c_{2}}{+}\delta)/2>0$. There is a $\rho_{2}>0$ such that {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}<\rho_{2} implies

[TABLE]

Now using the continuity of $\gamma(D_{1},D_{2})(\mu)$ as function of $\mu$ to combine these results, we find $\gamma(D_{1},D_{2})(\mu)$ has a zero in the interval $[\mu_{c_{2}}-\delta,\mu_{c_{2}}+\delta]$, provided that {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}<\min\{\rho_{1},\rho_{2}\}.

Combining the results of the previous two paragraphs, on the interval $[\mu_{c_{2}}-\delta,\mu_{c_{2}}+\delta]$ and for $(D_{1},D_{2})$ sufficiently close to $(D,D)$, the characteristic polynomial of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} has a complex conjugate pair of roots and at least one pair of purely imaginary roots. Moreover, by choice of $\delta$, $\gamma^{\prime}(D,D)(\mu)=A^{\prime}(\mu)>A^{\prime}(\mu_{c_{2}})/2>0$ .

To see that the transversality condition holds, let $\eta_{3}=A^{\prime}(\mu_{c_{2}})/2>0$. By lemma 3.7, if $(D_{1},D_{2})$ is sufficiently close to $(D,D)$, on the interval $[\mu_{c_{2}}-\delta,\mu_{c_{2}}+\delta]$ we have

[TABLE]

since $\gamma^{\prime}(D,D)(\mu)=A^{\prime}(\mu)>A^{\prime}(\mu_{c_{2}})/2$ on the interval $[\mu_{c_{2}}-\delta,\mu_{c_{2}}+\delta]$, in particular, at the point where $\gamma(D_{1},D_{2})(\mu)$ has a zero. This completes the proof that hypothesis 3 of theorem 3.3 holds for $E_{2}(\mu,D_{1},D_{2})$. ∎

4 Bifurcation to cycles

In [7] it is shown that cycles exist for certain values of parameters and under the assumption that $D{=}D_{1}{=}D_{2}$. In a nutshell, the assumption implies that the $\omega$-limit set of a solution starting near the unstable equilibrium $E_{2}(\mu,D,D)$ is contained in a plane in $NPZ$-space. Of course, this plane also contains $E_{2}$. The authors observe that there is a two-dimensional attracting set for this system. The Poincaré-Bendixson theorem applies to the two-dimensional limit system, delivering the existence of a cycle for the $(N,P,Z)$-system.

In this section we fix $(D_{1},D_{2})$ sufficiently close to $(D,D)$ so that theorem 3.8 applies, and we want to apply the Hopf bifurcation theorem [3, 8], restated in theorem 3.3, to deduce that the system (1) undergoes a Hopf bifurcation as $\mu$ passes a critical value. In section 3 we studied the characteristic polynomial of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} relative to the characteristic polynomial of J\bigl{(}E_{2}(\mu,D,D)\bigr{)}. Features of these polynomials continue to play a role. The results of section 3 show that for parameter values $(D_{1},D_{2})$ near $(D,D)$, hypotheses 3 and 4 of theorem 3.3 concerning the behavior of the eigenvalues of the Jacobian J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} are satisfied. However, in system (1), the coordinates of the equilibrium E_{2}\bigl{(}N(\mu,D_{1},D_{2}),\lambda_{Z}(D_{2}),Z(\mu,D_{1},D_{2})\bigr{)} are changing with the parameter $\mu$, so hypothesis 1 of theorem 3.3 is not satisfied. The immediate goal of this section is to overcome this difficulty by using the inverse function theorem.

Augment system (1) by introducing the bifurcation parameter as an extra variable, giving

[TABLE]

We are interested in the equilibria of system (23) as $\mu$ varies in a small interval around a number $\mu_{c_{2}}$ (depending on $D_{1}$ and $D_{2}$, but not necessarily uniquely determined) characterized in the proof of theorem 3.8 as a parameter value at which the Jacobian J\bigl{(}E(\mu_{c_{2}},D_{1},D_{2})\bigr{)} has a pair of purely imaginary eigenvalues, crossing the imaginary axis in the complex plane transversally. By lemma 3.1 the functions $N(\mu,D_{1},D_{2})$ and $Z(\mu,D_{1},D_{2})$ are smooth to the same degree of smoothness possessed by $f_{1}(N)$ and $f_{2}(P)$.

With $\delta$ as in the proof of theorem 3.8, let $I$ be an interval containing $\mu_{c_{2}}$, contained in $[\mu_{c_{2}}-\delta,\mu_{c_{2}}+\delta]$, and supporting a curve

[TABLE]

parameterizing the equilibrium locus of the augmented system (23) near the critical point $h(\mu_{c_{2}})$. We note that

[TABLE]

so $h$ is an immersion of an interval into ${\mathbf{R}}^{4}$.

Consider next the map $H\colon{\mathbf{R}}^{3}\times I\rightarrow{\mathbf{R}}^{4}$ defined by

[TABLE]

Observe that H(0,0,0,\mu)=\bigl{(}N(\mu,D_{1},D_{2}),\lambda_{Z}(D_{2}),Z(\mu,D_{1},D_{2}),\mu\bigr{)} and that

[TABLE]

which is an invertible matrix for any $\mu$. In particular, on a neighborhood $U$ of $(0,0,0,\mu)$, $H$ is a diffeomorphism of $U$ onto its image, smooth to the degree of smoothness of $N$ and $Z$, by the inverse function theorem [5, p.122]. Consequently, an interval of the form $(0,0,0){\times}I^{\prime}$ is mapped smoothly and bijectively onto the equilibrium locus of the system (23).

For $Y=(N,P,Z,\mu)$ in $H(U)$ we can write $Y=H(X)$, where $X=(x,y,z,\mu)$. To simplify the notation, write $dY/dt=F(Y)$ for system (23). Then

[TABLE]

is a formal expression for the system with respect to the new coordinates. If $Y_{0}$ is an equilibrium solution of system (23), and $H(X_{0})=Y_{0}$, then $X_{0}$ is an equilibrium solution of the system (24).

Let us now examine the relation between $DF(Y_{0})$, the Jacobian of system (23) at $Y_{0}$, and the Jacobian of the system (24) at $X_{0}$. We compute

[TABLE]

Thus, the Jacobian of system (24) at the equilibrium $X_{0}$ is simply a conjugate of the Jacobian of system (23) at the equilibrium $Y_{0}=H(X_{0})$. In particular, the eigenvalues of the Jacobians are the same. We can now prove our main theorem.

Theorem 4.1**.**

Assume the hypotheses of theorem 3.8 hold and that $f_{1}$ and $f_{2}$ are four times continuously differentiable. For $D_{1}$ and $D_{2}$ both sufficiently close to $D$ there is a value $\mu_{c_{2}}$ of the growth parameter at which the system (1) undergoes a Hopf bifurcation, resulting in the appearance of a cycle.

Proof.

Let us now make the assumption that {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)} is so small that the conclusions of theorem 3.8 hold, giving us a number $\mu_{c_{2}}$ (depending on $D_{1}$ and $D_{2}$, but not necessarily uniquely determined) at which the Jacobian J\bigl{(}E(\mu_{c_{2}},D_{1},D_{2})\bigr{)} has a purely imaginary pair of eigenvalues, crossing the imaginary axis in the complex plane transversally. The assumption that $f_{1}$ and $f_{2}$ are four times continuously differentiable fulfills part two of theorem 3.3. Now review elements of the proof pertaining to the eigenvalues of the Jacobian of system (23) at $(N(\mu,D_{1},D_{2}),\lambda_{Z}(D_{2}),Z(\mu,D_{1},D_{2}),\mu)$ as $\mu$ ranges over a small interval and increases through $\mu_{c_{2}}$. Throughout, the Jacobian has a negative real eigenvalue, by theorem 3.8, part one. By part two of theorem 3.8, for $\mu<\mu_{c_{2}}$ and sufficiently near $\mu_{c_{2}}$, there is a complex-conjugate pair of eigenvalues with negative real part. At $\mu=\mu_{c_{2}}$, the real part vanishes, and, for $\mu>\mu_{c_{2}}$ and sufficiently near $\mu_{c_{2}}$, there is a complex-conjugate pair of eigenvalues with positive real part. Moreover, the derivative of the function selecting the real part of the complex-conjugate pair is positive at $\mu_{c_{2}}$. That is, the locus of the complex-conjugate pair of eigenvalues crosses the imaginary axis transversally.

By our observations on the relationship of the system (23) to the system (24), as the parameter $\mu$ varies, the evolution of the eigenvalues of equilibria of system (24) has the same characteristics. Thus, hypotheses 4 and 3 of theorem 3.3 are satisfied for system (24). Having assumed that $f_{1}$ and $f_{2}$ are four times continuously differentiable, then hypothesis 2 of theorem 3.3 is also satisfied. Finally, since $(0,0,0)$ is an isolated equilibrium for all relevant values of $\mu$ for system (24), hypothesis 1 of theorem 3.3 is satisfied. The consequence is that the equilibrium $(0,0,0)$ undergoes a Hopf bifurcation at $\mu=\mu_{c_{2}}$, after which cycles appear in the phase portrait of system (24). Then the local diffeomorphism $H$ carries this portrait forward into the phase portrait of system (23). Thus, we have proved that cycles appear in our original system (1) when the parameter $\mu$ slightly exceeds $\mu_{c_{2}}$. ∎

One can make similar observations viewing the Routh-Hurwitz expressions $a_{1}$, $a_{2}$, and $a_{3}$ as functions of $(D_{1},D_{2})$, but it seems that this only delivers “loss of stability,” which is not quite enough. One can also play around with derivatives of J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)} with respect to $D_{1}$ and $D_{2}$. One might get information about how far $(D_{1},D_{2})$ can vary from $(D,D)$ and still have a result.

5 Simulations and Examples

We consider now system (1) and take $f_{1}$ and $f_{2}$ to be Michaelis-Menten functions, also called Holling type II functions. We choose notations as follows.

[TABLE]

With these choices explicit formulas can be given for many quantities studied in earlier sections. For example, we obtain formulas for the numbers $\lambda_{P}(D_{1})$ and $\lambda_{Z}(D_{2})$ defined in equation (2). For $\lambda_{P}(D_{1})$ we have

[TABLE]

for $\lambda_{Z}(D_{2})$ we have

[TABLE]

For rate functions of this type, the equation (6) determining $N(\mu,D_{1},D_{2})$ reduces to a quadratic equation. In principle, one obtains explicit solutions for $N(\mu,D_{1},D_{2})$ and the nonnegative solution is easily identified. Then the equation (7) for $Z(\mu,D_{1},D_{2})$ is linear and easily solved for $Z(\mu,D_{1},D_{2})$.

With these remarks we can turn to an illustration of theorems 3.5 and 3.8 in a case using Holling type II rate functions. With these rate functions, the concavity of the function $f_{2}$ guarantees that the condition D/\bigl{(}\gamma_{2}\lambda_{Z}(D)\bigr{)}>f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)} is satisfied. First, we set $D{=}D_{1}{=}D_{2}{=}1.0$ and the remaining parameters to the following values.

[TABLE]

With all quantities involved explicitly computed, the formula given in (20) for the real part of the complex conjugate pair of eigenvalues of the linearization of the system at the coexistence equilibrium can be made explicit, though messy, and is easily plotted by a computer algebra system. This is the solid curve in figure 2, which shows that a Hopf bifurcation occurs in the vicinity of $\mu=0.6$. Figure 3 exhibits solutions of the system for $\mu$ slightly smaller and slightly larger than the bifurcation value $\mu_{c_{2}}(D,D)\approx 0.6$.

Next, keep $D{=}1$ and set $D_{1}{=}1.2$ and $D_{2}{=}1.3$. Interpreted graphically, theorem 3.8 says that the graph of the function defined by taking the real part of the complex conjugate pair of eigenvalues associated with the coexistence equilibrium is an increasing function whose graph lies in a neighborhood of the curve we discussed in the preceding paragraph. We do not have an explicit formula for this function with the new values for $D_{1}$ and $D_{2}$, but we can can estimate its values by numerically computing the eigenvalues of the linearization along a sequence of $\mu$-values. Figure 2 also shows a sequence of points derived from eigenvalue approximations when $D_{1}{=}1.2$ and $D_{2}{=}1.3$. Interpolating a curve through the plotted points, we see a Hopf bifurcation occurs in the vicinity of $\mu=0.9$. Figure 4 shows trajectories of this system for values of $\mu$ slightly smaller and slightly larger than the bifurcation value.

To provide an additional illustration of the result of theorem 3.8, we consider a version of system (1) incorporating rate functions with the property that the graphs have inflection points. Consider

[TABLE]

with the parameter values

[TABLE]

Further, set

[TABLE]

Then the condition D/\bigl{(}\gamma_{2}\lambda_{Z}(D)\bigr{)}>f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)} of theorem 3.8 is satisfied, but this is not a consequence of the concavity of the graph of $f_{2}$.

Determining $\lambda_{P}(D_{1})$ and $\lambda_{Z}(D_{2})$ in this example requires solving quadratic equations, so obtaining exact values is quite easy. Consequently, one can compute explicitly from (5) the value $\mu_{c_{1}}(D_{1},D_{2})$ beyond which the coexistence equilibrium exists. Locating a coexistence equilibrium $E_{2}(\mu,D_{1},D_{2})$ requires solving a cubic equation for $N(\mu,D_{1},D_{2})$, for which it is more appropriate to use numerical methods. We choose a sequence of $\mu$ values starting beyond $\mu_{c_{1}}(D_{1},D_{2})$, approximate the coexistence equilibria and their Jacobian matrices J\bigl{(}E_{2}(\mu,D_{1},D_{2})\bigr{)}, and numerically compute the real part of the complex pair of eigenvalues of the Jacobian for each $\mu$ value. Plotting the real part against $\mu$ produces figure 5, which exhibits the expected change of sign and shows that a Hopf bifurcation occurs in the vicinity of $\mu=7.25$; trajectories for parameters slightly smaller and slightly larger than the bifurcation value are shown in figure 6.

6 Uniform Approximation

For the proof of theorem 3.8, we need lemma 3.7, which states that, when $(D_{1},D_{2})$ is close to $(D,D)$, $\gamma^{\prime}(D_{1},D_{2})(\mu)$ is uniformly approximated by $\gamma^{\prime}(D,D)(\mu)=A^{\prime}(\mu)$ on an interval $[\mu_{c_{2}}-\delta_{0},\mu_{c_{2}}+\delta_{0}]$, where $\mu_{c_{2}}$ is a point where $A(\mu_{c_{2}}){=}0$ and $A^{\prime}(\mu_{c_{2}})>0$.

To recapitulate theorem 3.5, the working assumptions are that a value $D$ is fixed, a coexistence equilibrium $E_{2}(\mu,D,D)$ exists, and that there is a range of parameters $\mu$ for which the linearizations at the coexistence equilibrium have a complex conjugate pair of eigenvalues. Moreover, at the parameter value $\mu_{c_{2}}$, the linearization of the system has a purely imaginary pair of eigenvalues. For any $\mu$ slightly smaller than $\mu_{c_{2}}$, the pair of complex eigenvalues has a negative real part and for any $\mu$ slightly larger than $\mu_{c_{2}}$ the pair of complex eigenvalues has a positive real part. By lemma 3.1 there is an interval $I$ containing $\mu_{c_{2}}$ and a disc $\Delta$ centered at $(D,D)$ such that the functions $N(\mu,D_{1},D_{2})$ and $Z(\mu,D_{1},D_{2})$ defined on $I{\times}\Delta$ smoothly parametrize the locus of coexistence equilibria near $E_{2}(\mu,D,D)$.

Uniform approximation of $\gamma^{\prime}(D_{1},D_{2})(\mu)$ by $\gamma^{\prime}(D,D)(\mu)$ implies that, for the coexistence equilibrium $E_{2}(\mu,D_{1},D_{2})$, there is a $\mu$-interval in $I$ for which the eigenvalues of the linearizations exhibit the same qualitative behavior as described for $E_{2}(\mu,D,D)$, as shown in theorem 3.8.

Proposition 6.1** (Lemma 3.7).**

Assume $f_{1}$ is three times continuously differentiable and $f_{2}$ is two times continuously differentiable. Then there exists a $\mu$-interval $[\mu{-}\hat{\delta},\mu{+}\hat{\delta}]$ on which the $\mu$-derivative $\gamma^{\prime}(D_{1},D_{2})(\mu)$ is uniformly approximated by $\gamma^{\prime}(D,D)(\mu)=A^{\prime}(\mu)$. In fact, there exists a constant $C$ such that

[TABLE]

for any $\mu\in[\mu_{c_{2}}{-}\hat{\delta},\mu_{c_{2}}{+}\hat{\delta}]$.

The proposition follows from a sequence of lemmas and estimates, given below. In the course of proving these results, we find it necessary to impose the differentiability conditions on $f_{1}$ and $f_{2}$.

An essential ingredient in the process is to obtain bounds on magnitudes of the differences

[TABLE]

in terms of {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}, and where the $p_{i}(\mu,D_{1},D_{2})=-a_{3-i}(\mu,D_{1},D_{2})$, $0\leq i\leq 2$ are given by the formulas (12), (13), and (14).

We now explain the role played by these bounds. By the chain rule, we compute $\gamma^{\prime}(D_{1},D_{2})(\mu)$ as the inner product of a row vector $\nabla\gamma$ with a column vector $(p_{0}^{\prime},p_{1}^{\prime},p_{2}^{\prime})$:

[TABLE]

Remark*.*

In order to avoid extremely long expressions in the following analysis, we use abbreviations such as

[TABLE]

We can write

[TABLE]

where $\langle 0,0,1\rangle$ represents the projection $\pi_{3}$ to the third coordinate $\gamma$ in $P_{1}(x)^{-}{\times}P_{2}(x)^{-}$. The gradients $\nabla\alpha$, $\nabla\beta$, and $\nabla\gamma$ are evaluated at

[TABLE]

where explicit expressions for $a_{i}(\mu,D_{1},D_{2})$ are given in (12), (13), and (14). The derivatives $p_{0}^{\prime}$, $p_{1}^{\prime}$, and $p_{2}^{\prime}$ are evaluated at $(\mu,D_{1},D_{2})$.

Applying the convention of (30), we have

[TABLE]

Before we go farther, we compress taking the derivative $D(M^{-1})$ at $(p_{0},p_{1},p_{2})$ and evaluating on the vector $(p_{0}^{\prime},p_{1}^{\prime},p_{2}^{\prime})$, writing

[TABLE]

Then the previous equation becomes

[TABLE]

We can estimate using operator norms computed in terms of the Euclidean metrics.

[TABLE]

since the norm of a projection is $1$. Now we use the triangle inequality to bound the last expression.

[TABLE]

Explicit expansion of the first summand in (32) is given in the proof of lemma 6.2, where we will see the role of the bounds on

[TABLE]

Similarly, explicit expansion of the second summand in (32) is given in the proof of lemma 6.3, where we will see the role of the bounds on

[TABLE]

Lemma 6.2**.**

There is a constant $C_{3}$ such that the first summand in the expansion (32) satisfies

[TABLE]

Proof.

By the basic property of the operator norm,

[TABLE]

Examining the factor coming from the operator norm,

[TABLE]

for some constant $C_{DM^{-1}}$. This is because, as $(D_{1},D_{2})$ ranges over any closed disc centered on $(D,D)$ and $\mu$ ranges over any closed interval containing $\mu_{c_{2}}$, the coordinates $(p_{0},p_{1},p_{2})$ are contained in a compact set, so there is a constant $C_{DM^{-1}}$ that bounds the norm of $D(M^{-1})$ at any of these points.

Now we require a bound on the other factor, for which we have

[TABLE]

provided by combining the results of propositions 6.5, 6.6, and 6.7. Combining this bound with the bound on the expression (34), the first summand in (32) is bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. ∎

Now we bound the second summand in (32).

Lemma 6.3**.**

There is a constant $C_{4}$ such that the second summand in the expansion (32) satisfies

[TABLE]

Proof.

To handle the second summand in (32), we have the bound

[TABLE]

To the first factor in the bounding term, apply the mean value theorem [5, p.103, Corollary 1] for vector-valued functions of several variables, obtaining

[TABLE]

The derivative of the map $P_{3}(x)^{-}\rightarrow M_{3,3}({\mathbf{R}})$, $(p_{0},p_{1},p_{2})\mapsto D(M^{-1})(p_{0},p_{1},p_{2})$ is continuous, because $M^{-1}$ is $C^{\infty}$. The evaluation point $(q_{0},q_{1},q_{2})$ is on the line segment connecting $[(p_{0},p_{1},p_{2})](\mu,D_{1},D_{2})$ and $[(p_{0},p_{1},p_{2})](\mu,D,D)$. Again the possibilities range over a compact set, so the norm of the second derivative satifies

[TABLE]

for a constant $C_{D^{2}(M^{-1})}$ independent of $(D_{1},D_{2},\mu)$. We compute

[TABLE]

by combining the results of propositions 6.5, 6.6, and 6.7. Assembling the bounds in (38) and (39), \bigl{\lVert}[D(M^{-1})](\mu,D_{1},D_{2})-[D(M^{-1})](D,D,\mu)\bigr{\rVert} is bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. This takes care of the first factor on the righthand side of (37).

For the remaining factor in the bounding term in (37), the norm of the tangent vector satisfies

[TABLE]

for some constant $C_{T}$, independent of $\mu$. Now that we have taken care of both factors on the righthand side of (37) the second summand in (32) is bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}, as claimed. ∎

We can now prove the uniform convergence result.

Proof of proposition 6.1..

Combining the inequalities (31) and (32), we have

[TABLE]

Using the observations detailed in lemmas 6.2 and 6.3,

[TABLE]

Thus, there is a constant $C$ such that

[TABLE]

for any $\mu\in[\mu_{c_{2}}-\delta_{0},\mu_{c_{2}}+\delta_{0}]$. ∎

We have already used a technique of obtaining bounds by splitting quantities. As we will continue to exploit the technique in the following results, we formulate the lemma 6.4 for reference.

Lemma 6.4**.**

Let $Q_{1}(\mu,D_{1},D_{2})$ and $Q_{2}(\mu,D_{1},D_{2})$ be quantities defined on a domain $I{\times}\Delta$ satisfying the following conditions.

There is a constant $c_{1}$ such that

[TABLE] 2. 2.

There is a constant $c_{2}$ such that

[TABLE] 3. 3.

There are constants $c_{3}$ and $c_{4}$ such that $\lvert Q_{1}(\mu,D_{1},D_{2})\rvert\leq c_{3}$ for $(\mu,D_{1},D_{2})\in I{\times}\Delta$ and for {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)} sufficiently small, and $\lvert Q_{2}(\mu,D,D)\rvert\leq c_{4}$ for $\mu\in I$.

Then there is a constant $c_{5}$ such that, for {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)} sufficiently small,

[TABLE]

As has been seen, the proof of proposition 6.1 depends on the following three propositions.

Proposition 6.5**.**

There are constants $C_{0}$ and $C_{0}^{\prime}$ such that

[TABLE]

Proposition 6.6**.**

There are constants $C_{1}$ and $C_{1}^{\prime}$ such that

[TABLE]

Proposition 6.7**.**

There are constants $C_{2}$ and $C_{2}^{\prime}$ such that

[TABLE]

The proofs of propositions 6.5, 6.6, and 6.7 depend in turn on a number of elementary bounds and estimates, given below in lemma 6.8 and propositions 6.9, 6.10, and 6.11. We give the quick proofs of lemma 6.8 and propositions 6.9 and 6.10, because they are quite short, postponing the proof of the many parts of proposition 6.11 to the end of the section. After we state these results, we prove propositions 6.5, 6.6, and 6.7.

Lemma 6.8**.**

Given $D$, there is an interval $J$ containing $\lambda_{Z}(D)$ such that

[TABLE]

for all $P\in J$. Thus, for all $P\in J$, $f_{2}^{\prime}(P)$ is bounded away from zero, and, for all $D_{2}$ in the preimage $\lambda_{Z}^{-1}(J)$, f_{2}^{\prime}\bigl{(}\lambda_{Z}(D_{2})\bigr{)} is bounded away from zero.

Proof.

By assumption on $f_{2}$, f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)}>0. By continuity of $f_{2}^{\prime}$, there is an interval $J$ containing $\lambda_{Z}(D)$ such that, for all $P\in J$,

[TABLE]

Proposition 6.9**.**

Each of the quantities

[TABLE]

is bounded by some constant on the interval $I$.

Proof.

Each of the listed functions is continuous on the closed interval $I$, so each one is bounded. ∎

Proposition 6.10**.**

Each of the quantities

[TABLE]

is bounded by some constant on the domain $I{\times}\Delta$.

Proof.

Each of the listed functions is continuous on the compact set $I{\times}\Delta$, so each one is bounded. ∎

The proof of the next result depends on many more details of system (1) and consequences drawn from them. Some steps in the proof are quite lengthy, so we postpone these details to the end of the section.

Proposition 6.11**.**

Assume that the domain $I{\times}\Delta$ for which $N(\mu,D_{1},D_{2})$ and $Z(\mu,D_{1},D_{2})$ are defined is a subset of $I{\times}{\mathbf{R}}{\times}\lambda_{Z}^{-1}(J)$, where $J$ is as in lemma 6.8. Then for $\mu\in I$ and $(D_{1},D_{2})\in\Delta$, these differences are bounded by constants times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}.

[TABLE]

where the derivatives are taken with respect to $\mu$, are also bounded by constants times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}.

Proof of proposition 6.5.

After some reorganization, we have from (14)

[TABLE]

To bound $\lvert p_{0}(\mu,D_{1},D_{2})-p_{0}(D,D,\mu)\rvert$, we bound the absolute values of summands in (47) as follows. First, note \lvert D{-}D_{2}\rvert\leq{\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}, so we bound

[TABLE]

by using lemma 6.4, lemma 6.9, and lemma 6.10 to combine the noted bound with bounds 4 and 5 from proposition 6.11; bound

[TABLE]

using lemma 6.4 to combine bounds 1, 3, 4 and 5 with \lvert D{-}D_{2}\rvert\leq{\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. We compute from (47)

[TABLE]

For $\lvert p_{0}^{\prime}(\mu,D_{1},D_{2})-p_{0}^{\prime}(D,D,\mu)\rvert$, we bound from the first line of the expansion

[TABLE]

by combining bounds 5 and 8 with the bound \lvert D{-}D_{2}\rvert\leq{\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}; bound from the second and third lines

[TABLE]

by combining bounds 1, 3, 5, and 8 with the bound \lvert D{-}D_{2}\rvert\leq{\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}; bound from the fourth and fifth lines

[TABLE]

by combining bounds 1, 4, 5, 6, and 7 from proposition 6.11 with the bound \lvert D{-}D_{2}\rvert\leq{\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. ∎

Proof of proposition 6.6.

After some reorganization, placing terms belonging to $p_{1}(\mu,D_{1},D_{2})$ down the left side of the display, we have from (13)

[TABLE]

To bound $\lvert p_{1}(\mu,D_{1},D_{2})-p_{1}(D,D,\mu)\rvert$ by a constant multiple of {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}, first observe that both $\lvert D{-}D_{1}\rvert$ and $\lvert D{-}D_{2}\rvert$ are bounded by {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. Similarly bound

[TABLE]

by combining bounds 1, 4, 3, and 5 from proposition 6.11; bound

[TABLE]

is taken care of in bound 2 of proposition 6.11. Finally, $\lvert-DD_{1}+D^{2}\rvert$ has already been taken care of. Adding all these bounds, $\lvert p_{1}(\mu,D_{1},D_{2})-p_{1}(D,D,\mu)\rvert$ is bounded by a constant multiple of {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}.

Rather than exhibit a complete formula for $p_{1}^{\prime}(\mu,D_{1},D_{2})-p_{1}^{\prime}(D,D,\mu)$, we pick apart equation (48) to express this difference as a sum of expressions. From the first two lines,

[TABLE]

is involved in the sum. From lines three through six, the expressions involved are

[TABLE]

Making several applications of lemma 6.4, lemma 6.9, lemma 6.10, and proposition 6.11, we find that the absolute value of each of the quantities displayed in (49), (50), (51), (52), and (53) is bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. Consequently, there is a constant $C_{1}^{\prime}$ such that

[TABLE]

Proof of proposition 6.7.

Using the formula (12) and organizing the difference $p_{2}(\mu,D_{1},D_{2})-p_{2}(\mu,D,D)$ to display terms belonging to $p_{2}(\mu,D_{1},D_{2})$ down the left side of the display, we have

[TABLE]

The bound on $\lvert p_{2}(\mu,D_{1},D_{2})-p_{2}(\mu,D,D)\rvert$ in (45) follows from bounds on summands in (54), as follows. Apply lemma 6.4, lemma 6.9, lemma 6.10, and bounds 4 and 5 from proposition 6.11 to bound the term

[TABLE]

using bound 2 from proposition 6.11. Finally, $\lvert-D_{1}+D\rvert$ is bounded by {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}.

From (54),

[TABLE]

For (46), the same sort of building block approach on terms in (55) using bounds 5 and 8 from proposition 6.11, then bounds 1, 6, and 7, and finally bounds 3 and 7 delivers the bound (46). ∎

Now we embark on the proof of proposition 6.11. Part of this work is made easier by the fact that the rate functions in (1) are explicitly linear in $D$, $D_{1}$, and $D_{2}$. On the other hand, because other quantities such as $N$, $P$, and $Z$ depend implicitly on $D$, $D_{1}$, and $D_{2}$, the details of the analyses are somewhat lengthy.

Proof of bound 1.

By definition and the mean value theorem applied to $f_{2}$, we have

[TABLE]

for some number $P_{1}$ between $\lambda_{Z}(D_{2})$ and $\lambda_{Z}(D)$. Consequently,

[TABLE]

Since $P_{1}$ is also close to $\lambda_{Z}(D)$, we may assume by lemma 6.8 that f_{2}^{\prime}(P_{1})>f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)}/2. Thus, for $D_{2}$ sufficiently close to $D$

[TABLE]

Proof of bound 2.

Fix $\mu$ in the interval $I$. We may use the mean-value theorem [5, p.103, Corollary 1] for functions of variables $(D_{1},D_{2})$, obtaining

[TABLE]

where $(\widehat{D}_{1},\widehat{D}_{2})$ is a point on the line segment connecting $(D_{1},D_{2})$ and $(D,D)$. Therefore, we have to bound the magnitude of the gradient $\lVert\nabla(f_{1}\circ N)(\widehat{D}_{1},\widehat{D}_{2})\rVert$ in a disc surrounding $(D,D)$ by a constant.

To obtain information about the partial derivatives $\partial(f_{1}\circ N)/\partial D_{1}=f_{1}^{\prime}(N)\cdot(\partial N/\partial D_{1})$ and $\partial(f_{1}\circ N)/\partial D_{2}=f_{1}^{\prime}(N)\cdot(\partial N/\partial D_{2})$, we return to the defining equation

[TABLE]

and differentiate with respect to $D_{1}$ and $D_{2}$. We obtain

[TABLE]

From the first of these equations

[TABLE]

Now we obtain a bound on the gradient via

[TABLE]

since f^{\prime}_{2}\bigl{(}\lambda_{Z}(D_{2})\bigr{)}>f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)}/2 by choice of $\Delta$ and lemma 6.8. Also, f_{1}\bigl{(}(N(\mu,D_{1},D_{2})\bigr{)} is bounded by $\lim_{N\rightarrow\infty}f_{1}(N)$, since $N(\mu,D_{1},D_{2})$ is unbounded as $\mu$ tends to infinity according to theorem 3.2. Thus, $\lvert\partial N/\partial D_{2}(\mu,D_{1},D_{2})\rvert$ is bounded by a constant in the disc $\Delta$.

We also observe that f_{1}^{\prime}\bigl{(}N(\mu,\widehat{D}_{1},\widehat{D}_{2})\bigr{)} is bounded by a constant depending only on $I{\times}\Delta$, because of the convexity of the closed disc $\Delta$ centered at $(D,D)$ from which we choose $(D_{1},D_{2})$.

Combining all this information, $\lVert\nabla(f_{1}\circ N)(\widehat{D}_{1},\widehat{D}_{2})\rVert$ is bounded by a constant depending on $I{\times}\Delta$. Therefore, \bigl{\lvert}f_{1}\bigl{(}N(\mu,D_{1},D_{2})\bigr{)}-f_{1}\bigl{(}N(\mu,D,D)\bigr{)}\bigr{\rvert} is bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. ∎

Proof of bound 3.

We again use the mean-value theorem [5, p.103, Corollary 1] for functions of variables $(D_{1},D_{2})$, obtaining

[TABLE]

where $(\widehat{D}_{1},\widehat{D}_{2})$ is a point on the line segment connecting $(D_{1},D_{2})$ and $(D,D)$. Therefore, we have to bound the magnitude of the gradient $\lVert\nabla(f^{\prime}_{1}\circ N)(\widehat{D}_{1},\widehat{D}_{2})\rVert$ in a disc surrounding $(D,D)$.

To obtain information about the partial derivatives $\partial(f_{1}^{\prime}\circ N)/\partial D_{1}=f_{1}^{(2)}(N)\cdot(\partial N/\partial D_{1})$ and $\partial(f_{1}^{\prime}\circ N)/\partial D_{2}=f_{1}^{(2)}(N)\cdot(\partial N/\partial D_{2})$, we cite (56) for the vanishing of $\partial N/\partial D_{1}$ and the bound on $\partial N/\partial D_{2}$ obtained in (59). We also observe that f_{1}^{(2)}\bigl{(}N(\mu,\widehat{D}_{1},\widehat{D}_{2})\bigr{)} is bounded by a constant, assuming a continuous second derivative of $f_{1}$. The constant depends on the $\mu$-interval $I$ and the closed disc $\Delta$ centered at $(D,D)$ from which we choose $(D_{1},D_{2})$, plus the convexity of the disc.

Combining all this information, $\lVert\nabla(f_{1}^{\prime}\circ N)(\widehat{D}_{1},\widehat{D}_{2})\rVert$ is bounded by a constant depending on $I{\times}\Delta$. Therefore, \bigl{\lvert}f_{1}^{\prime}\bigl{(}N(\mu,D_{1},D_{2})\bigr{)}-f_{1}^{\prime}\bigl{(}N(\mu,D,D)\bigr{)}\bigr{\rvert} is bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. ∎

Proof of bound 4.

Since $\mu$ is a fixed number in $I$, we use again the mean value theorem for functions of $(D_{1},D_{2})$.

[TABLE]

where $(\widehat{D}_{1},\widehat{D}_{2})$ is a point on the line segment connecting $(D_{1},D_{2})$ and $(D,D)$. Therefore, we have to bound the magnitude of the gradient $\lVert\nabla Z(\widehat{D}_{1},\widehat{D}_{2})\rVert$ in a disc surrounding $(D,D)$.

To obtain information about $\partial Z/\partial D_{1}$ and $\partial Z/\partial D_{2}$ we return to the defining equation

[TABLE]

Differentiating with respect to $D_{1}$, we obtain

[TABLE]

since $\partial N/\partial D_{1}=0$ by equation (56). Differentiating with respect to $D_{2}$, we obtain

[TABLE]

Rewriting these equations, we obtain

[TABLE]

To bound $\partial Z/\partial D_{1}$, we require a bound on $\lambda_{Z}(D_{2})$. By definition f_{2}\bigl{(}\lambda_{Z}(D_{2})\bigr{)}=D_{2}/\gamma_{2}, so restricting $D_{2}$ to be close to $D$ prevents $D_{2}/\gamma_{2}$ from approaching $\lim_{P\rightarrow\infty}f_{2}(P)$. Consequently, $\lambda_{Z}(D_{2})$ is a bounded distance from $\lambda_{Z}(D)$.

To bound $\partial Z/\partial D_{2}$, we discuss the terms on the righthand side of (61) in reverse order. The term $D_{2}^{-1}Z(\mu,D_{1},D_{2})$ is bounded if $(D_{1},D_{2})$ is close to $(D,D)$, for it will be close to $D_{2}^{-1}Z(\mu,D,D)$. In turn $Z(\mu,D,D)$ is bounded by $\lim_{\mu\rightarrow\infty}Z(\mu,D,D)$, which exists by theorem 3.2. Concerning the second term, f_{2}^{\prime}\bigl{(}\lambda_{Z}(D_{2})\bigr{)}\cdot\lambda_{Z}^{\prime}(D_{2})=1/\gamma_{2}, so

[TABLE]

where f_{1}\bigl{(}N(\mu,D_{1},D_{2})\bigr{)} can be bounded in terms of $\lim_{N\rightarrow\infty}f_{1}(N)$, and f_{2}^{\prime}\bigl{(}\lambda_{Z}(D_{2})\bigr{)} is bounded away from zero by lemma 6.8. Concerning the first term, substitute the expression for $\partial N/\partial D_{2}$ given in (57), obtaining

[TABLE]

Consequently,

[TABLE]

where we use again the fact that \lambda_{Z}^{\prime}(D_{2})=1/\gamma_{2}f_{2}^{\prime}\bigl{(}\lambda_{Z}(D_{2})\bigr{)}. Arguing as above, we conclude this term can be bounded by a constant, and, therefore, $\lvert\partial Z/\partial D_{2}(\mu,D_{1},D_{2})\rvert$ itself is bounded by a constant depending only on the domain $I{\times}\Delta$.

Combining these bounds $\lVert\nabla Z(\widehat{D}_{1},\widehat{D}_{2})\rVert$ is bounded by a constant, so we conclude that $\lvert Z(\mu,D_{1},D_{2})-Z(\mu,D,D)\rvert$ is bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. ∎

Proof of bound 5.

By the mean value theorem

[TABLE]

for some $P_{2}$ between $\lambda_{Z}(D_{2})$ and $\lambda_{Z}(D)$. Moreover,

[TABLE]

for some $P_{1}^{\prime}$ between $\lambda_{Z}(D_{2})$ and $\lambda_{Z}(D)$. We may combine to obtain

[TABLE]

Since we have control of the continuous derivatives $f_{2}^{\prime}$ and $f_{2}^{(2)}$ on the interval $J$ around $\lambda_{Z}(D)$, the difference \lvert f_{2}^{\prime}\bigl{(}\lambda_{Z}(D_{2})\bigr{)}-f_{2}^{\prime}\bigl{(}\lambda_{Z}(D)\bigr{)}\rvert is indeed bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. ∎

Proof of bound 6.

This is precisely parallel to the proofs of bounds 2 and 3. We may use the mean-value theorem for functions of variables $(D_{1},D_{2})$, obtaining

[TABLE]

where $(\widehat{D}_{1},\widehat{D}_{2})$ is a point on the line segment connecting $(D_{1},D_{2})$ and $(D,D)$. Therefore, we have to bound the magnitude of the gradient $\lVert\nabla(f_{1}^{(2)}\circ N)(\widehat{D}_{1},\widehat{D}_{2})\rVert$ in a disc surrounding $(D,D)$.

To bound the partial derivatives

[TABLE]

we have the vanishing of $\partial N/\partial D_{1}$ by (56) and a bound on $\lvert\partial N/\partial D_{2}\rvert$ from (59). We also observe that f_{1}^{(3)}\bigl{(}N(\mu,\widehat{D}_{1},\widehat{D}_{2})\bigr{)} is bounded by a constant, assuming a continuous third derivative of $f_{1}$.

Combining all this information, $\lVert\nabla(f_{1}^{(2)}\circ N)(\widehat{D}_{1},\widehat{D}_{2})\rVert$ is bounded by a constant depending on $I{\times}\Delta$. Therefore, \lvert f_{1}^{(2)}\bigl{(}N(\mu,D_{1},D_{2})\bigr{)}-f_{1}^{(2)}\bigl{(}N(\mu,D,D)\bigr{)}\rvert is bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. ∎

Proof of bound 7.

For this proof, return to the defining relation for $N(\mu,D_{1},D_{2})$, namely,

[TABLE]

and differentiate with respect to $\mu$, obtaining

[TABLE]

so

[TABLE]

Consequently,

[TABLE]

Applying lemma 6.4 to (62) with bounds 1 and 3 as input, we find that $\lvert N^{\prime}(\mu,D_{1},D_{2})-N^{\prime}(\mu,D,D)\rvert$ is bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. ∎

Proof of bound 8.

For this proof, return to the defining relation for $Z$, namely,

[TABLE]

and differentiate with respect to $\mu$, obtaining

[TABLE]

Thus,

[TABLE]

and

[TABLE]

Obviously \lvert D-D_{2}\rvert\leq{\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}, so we make several applications of lemma 6.4 to combine this fact with bounds 1, 3, and 7 to deduce that $\lvert Z^{\prime}(\mu,D_{1},D_{2})-Z^{\prime}(\mu,D,D)\rvert$ is bounded by a constant times {\rm dist}\bigl{(}(D_{1},D_{2}),(D,D)\bigr{)}. ∎