Phytoplankton Competition for Nutrients and Light in a Stratified Lake: A Mathematical Model Connecting Epilimnion and Hypolimnion

Abstract

A mathematical model connecting epilimnion and hypolimnion is proposed to describe the competition of phytoplankton for nutrients and light in a stratified lake. The existence and stability of nonnegative steady-state solutions are completely characterized for all possible parameter ranges by means of stability analysis, bifurcation theory, and extensive simulations. The critical thresholds for settling speed of phytoplankton cells in the thermocline and the loss rate of phytoplankton are established, which determine the survival or extirpation of phytoplankton in epilimnion and hypolimnion. In particular, it is shown that in two extreme cases, the principle of competitive exclusion always holds in a stratified lake. We also consider the influence of environmental parameters on the vertical distribution and biomass density of phytoplankton via a systematic sensitivity analysis, and investigate their roles in phytoplankton blooms. These results can be used for the prediction of phytoplankton competition and blooms in a stratified lake.

Appendices

Appendix A

In order to obtain the existence and stability of semi-trivial steady-state solutions of (2.1), we consider an eigenvalue problem

$$\begin{aligned} {\left\{ \begin{array}{ll} D\phi ''(x)=\lambda \phi (x),~x\in (0,x_h),\\ D\phi '(0)=\displaystyle \frac{\lambda }{\delta _1\lambda +\delta _2}\phi (0),~\phi (x_h)=0, \end{array}\right. } \end{aligned}$$

(A.1)

where $D,\delta _1,\delta _2>0$. A direct calculation shows that if $\lambda$ is a complex eigenvalue of (A.1), then ${\bar{\lambda }}$ is also a complex eigenvalue of (A.1).

Lemma A.1

If $\lambda$ is an eigenvalue of (A.1), then ${{\,\mathrm{Re}\,}}\lambda <0$.

Proof

If $\lambda$ is an eigenvalue (real or complex valued) and $\phi (x)$ is the corresponding eigenfunction, then

$$\begin{aligned} \lambda \displaystyle \int _0^{x_h}|\phi (x)|^2\mathrm{d}x= D\displaystyle \int _0^{x_h}\phi ''(x){\bar{\phi }}(x)\mathrm{d}x =-\displaystyle \frac{\lambda |\phi (0)|^2}{\delta _1\lambda +\delta _2}-D\displaystyle \int _0^{x_h} |\phi '(x)|^2\mathrm{d}x. \end{aligned}$$

Let

$$\begin{aligned} c_1=\displaystyle \int _0^{x_h}|\phi (x)|^2\mathrm{d}x,~c_2=D\displaystyle \int _0^{x_h} |\phi '(x)|^2\mathrm{d}x. \end{aligned}$$

Then, $\lambda$ satisfies

$$\begin{aligned} \delta _1c_1\lambda ^2+(\delta _2c_1+|\phi (0)|^2+\delta _1c_2)\lambda +\delta _2c_2=0. \end{aligned}$$

This means that ${{\,\mathrm{Re}\,}}\lambda <0$. This completes the proof. $\square$

Proof of Theorem 3.1

From (3.1), it is easy to see that $E_1\equiv (0,M_b,0,M_b)$ is the unique nutrients-only semi-trivial steady state solution of (2.1). The stability of $E_1$ is determined by the eigenvalue problem

$$\begin{aligned}&\lambda \xi =\left( \delta _a^*-\delta \right) \xi , \end{aligned}$$

(A.2a)

$$\begin{aligned}&\lambda \zeta =\theta \left( p\delta -\delta _a^*-\displaystyle \frac{a}{x_e}\right) \xi -\displaystyle \frac{b}{x_e}\zeta +\displaystyle \frac{b}{x_e}\psi (0),\end{aligned}$$

(A.2b)

$$\begin{aligned}&\lambda \varphi =D_b\varphi ''(x)-v\varphi '(x)+ \left( rf(M_b)g(I(x,0,0)) -\delta \right) \varphi ,~0<x<x_h,\end{aligned}$$

(A.2c)

$$\begin{aligned}&\lambda \psi =D_m\psi ''(x)+ \left( \theta p\delta -\theta rf(M_b)g(I(x,0,0))\right) \varphi ,~0<x<x_h,\end{aligned}$$

(A.2d)

$$\begin{aligned}&D_b\varphi '(0)-v\varphi (0)=-a\xi ,~D_b\varphi '(x_h)-v\varphi (x_h)=0, \end{aligned}$$

(A.2e)

$$\begin{aligned}&D_m\psi '(0)=b(\psi (0)-\zeta ),~\psi (x_h)=0. \end{aligned}$$

(A.2f)

To establish the local stability of $E_1$, we let $\lambda _1$ be the eigenvalue of (A.2) with the largest real part and let $(\xi ,\zeta ,\varphi ,\psi )$ be the corresponding eigenfunction. Note that the linearized system (A.2) is partially decoupled. We consider the following three cases: (i) $\xi \ne 0$; (ii) $\xi =0, \; \varphi \not \equiv 0$; or (iii) $\xi =0, \varphi \equiv 0$.

• Case (i): $\xi \ne 0$. In this case, the stability of $E_1$ is determined by (A.2a). Then, $\lambda _1=\delta _a^*-\delta$.

• Case (ii): $\xi =0, \; \varphi \not \equiv 0$. This case means that the stability of $E_1$ is determined by (A.2c) and its boundary condition (A.2e). Then, $\lambda _1=\delta _*-\delta$.

• Case (iii): $\xi =0, \; \varphi \equiv 0$. In this case, (A.2) reduces to

  $$\begin{aligned} {\left\{ \begin{array}{ll} \lambda \zeta =-\displaystyle \frac{b}{x_e}\zeta +\displaystyle \frac{b}{x_e}\psi (0),&{}\\ \lambda \psi (x)=D_m\psi ''(x),~~0<x<x_h,&{}\\ D_m\psi '(0)=b(\psi (0)-\zeta ),~\psi (x_h)=0.&{} \end{array}\right. } \end{aligned}$$

  (A.3)

  From the boundary condition of (A.3), we have $\zeta =0$ if $\psi \equiv 0$. On the other hand, if $\zeta =0$ in (A.3), we also have $\psi \equiv 0$. This shows that $\zeta \ne 0$ and $\psi (x)\not \equiv 0$. Then, either $\lambda =-b/x_e=-D_mk^2\pi ^2/x_h^2$ for some $k\in {\mathbb {N}}$ with the corresponding eigenfunction $(\zeta ,\psi )=(1,-\sqrt{bx_e/D_m}\sin (\sqrt{b/(x_eD_m)}x)$; or $\lambda \ne -b/x_e$, $\displaystyle \zeta =\displaystyle \frac{b}{x_e\lambda +b}\psi (0)$, and $\psi$ satisfies

  $$\begin{aligned} {\left\{ \begin{array}{ll} \lambda \psi (x)=D_m\psi ''(x),~~0<x<x_h,&{}\\ D_m\psi '(0)=\displaystyle \frac{bx_e\lambda }{x_e\lambda +b}\psi (0),~\psi (x_h)=0.&{} \end{array}\right. } \end{aligned}$$

  (A.4)

It follows from Lemma A.1 that all eigenvalues of (A.4) have negative real parts. Therefore, we conclude that in case (iii), ${{\,\mathrm{Re}\,}}\lambda _1$ is negative.

Summarizing the above cases (i)-(iii), we conclude that if (3.6) holds, then ${{\,\mathrm{Re}\,}}\lambda _1<0$ and $E_1$ is locally asymptotically stable; Conversely, if (3.7) holds, then ${{\,\mathrm{Re}\,}}\lambda _1>0$ and $E_1$ is unstable.

Next we prove that if (3.8) holds, then $E_1$ is globally asymptotically stable. From the first equation of (2.1), we have

$$\begin{aligned} \displaystyle \frac{dA(t)}{dt}\le rA\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,0))\mathrm{d}x -\displaystyle \frac{a}{x_e}A-\delta A, \end{aligned}$$

which implies that $\displaystyle \lim _{t\rightarrow \infty }A(t)=0$ if $\delta >\delta _a^{**}$. From the theory of asymptotical autonomous systems (Mischaikow et al. 1995), the third equation of (2.1) reduces to a limiting system

$$\begin{aligned} \displaystyle \frac{\partial B(x,t)}{\partial t}&=D_b\displaystyle \frac{\partial ^2B}{\partial x^2}-v\displaystyle \frac{\partial B}{\partial x}+rBf(M)g(I(x,A,B))-\delta B\\&\le D_b\displaystyle \frac{\partial ^2B}{\partial x^2}-v\displaystyle \frac{\partial B}{\partial x}+rg(I(0,0,0))B-\delta B,~0<x<x_h,~t>0,\\ D_b\displaystyle \frac{\partial B(0,t)}{\partial x}&-vB(0,t)= D_b\displaystyle \frac{\partial B(x_h,t)}{\partial x}-vB(x_h,t)=0,~t>0. \end{aligned}$$

This implies that B(x, t) converges to 0 uniformly for $x\in [0,x_h]$ as $t\rightarrow \infty$ if $\delta >\delta _{**}$ by the comparison principle of parabolic equations. By applying the theory of asymptotical autonomous systems again, (2.1) reduces to a limiting system

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{dN(t)}{dt}=\displaystyle \frac{b}{x_e}(M(0,t)-N(t)),&{}t>0,\\ \displaystyle \frac{\partial M(x,t)}{\partial t}=D_m\displaystyle \frac{\partial ^2M}{\partial x^2},&{}0<x<x_h,~t>0,\\ D_m\displaystyle \frac{\partial M(0,t)}{\partial x}=b(M(0,t)-N(t)),~M(x_h,t)=M_b,&{}t>0. \end{array}\right. } \end{aligned}$$

(A.5)

Define Lyapunov functional $V:{\mathbb {R}}\times C([0,L])\rightarrow {\mathbb {R}}$ by

$$\begin{aligned} V(N,M)=\displaystyle \frac{x_e}{2}(N-M_b)^2+\displaystyle \frac{1}{2}\displaystyle \int _0^{x_h}(M(z)-M_b)^2\mathrm{d}z. \end{aligned}$$

For an arbitrary solution (N(t), M(x, t)) of (A.5) with nonnegative initial values, we obtain

$$\begin{aligned} \displaystyle \frac{dV(N(t),M(x,t))}{dt}&=x_e(N(t)-M_b)\displaystyle \frac{d N}{dt} +\displaystyle \int _0^{x_h}(M(z,t)-M_b)\displaystyle \frac{\partial M}{\partial t}\mathrm{d}z\\&=b(N(t)-M_b)(M(0,t)-N(t))+D_m\displaystyle \int _0^{x_h}(M(z,t)-M_b)\displaystyle \frac{\partial ^2 M}{\partial z^2}\mathrm{d}z\\&=b(N(t)-M_b)(M(0,t)-N(t))\\&\quad +D_m\displaystyle \frac{\partial M}{\partial z}(M(z,t)-M_b)\Big |_0^{x_h}-\displaystyle \int _0^{x_h} \left( \displaystyle \frac{\partial M}{\partial z}\right) ^2\mathrm{d}z\\&=b(N(t)-M_b)(M(0,t)-N(t))\\&\quad -b(M(0,t)-N(t))(M(0,t)-M_b)-\displaystyle \int _0^{x_h} \left( \displaystyle \frac{\partial M}{\partial z}\right) ^2\mathrm{d}z\\&=-b(N(t)-M(0,t))^2-\displaystyle \int _0^{x_h} \left( \displaystyle \frac{\partial M}{\partial z}\right) ^2\mathrm{d}z\le 0. \end{aligned}$$

Note that $dV(\cdot )/dt=0$ holds if and only if $\partial M/\partial z=0$ and $N(t)=M(0,t)$. It follows from $M(x_h,t)=M_b$ that $N(t)\equiv M(x,t)\equiv M_b$. By using the LaSalle’s Invariance Principle (Henry 1981), we conclude that (N(t), M(x, t)) converges to $(M_b,M_b)$ uniformly for $x\in [0,x_h]$ as $t\rightarrow \infty$, and $E_1$ is globally asymptotically stable for (2.1) with respect to any nonnegative initial value. $\square$

Proof of Theorem 3.3

The steady-state equation (3.2) can be explicitly solved. The equation of N implies that $b(M(0)-N)=(1-p)\theta \delta x_eA$. Combining the equation of M with its boundary conditions $D_mM'(0)=(1-p)\theta \delta x_eA$ and $M(x_h)=M_b$ gives

$$\begin{aligned} M(x)=M_b-\displaystyle \frac{x_h-x}{D_m}(1-p)\theta \delta x_eA, \end{aligned}$$

and then

$$\begin{aligned} N=M_b-\left( \displaystyle \frac{x_h}{D_m}+\displaystyle \frac{1}{b}\right) (1-p)\theta \delta x_eA. \end{aligned}$$

If $p=1$ and (3.9) holds, then $M(x)=N=M_b$, and there exists a unique positive $A_2^*$ satisfying $\displaystyle f(M_b)\displaystyle \frac{r}{x_e}\int _{-x_e}^0g(I(x,A_2^*))\mathrm{d}x=\delta$ since g is a strictly decreasing function of A. If $p\in [0,1)$, then we let

$$\begin{aligned} Q=\delta A,~~\beta =\left( \displaystyle \frac{x_h}{D_m}+\displaystyle \frac{1}{b}\right) (1-p)\theta x_e. \end{aligned}$$

From the first equation of (3.2), we consider

$$\begin{aligned} u(\delta ,Q)= & {} \displaystyle \frac{r}{x_e}\displaystyle \int _{-x_e}^0g\left( I\left( x,\displaystyle \frac{Q}{\delta }\right) \right) \mathrm{d}x -\delta \left( 1+\displaystyle \frac{\gamma _1}{M_b-\beta Q}\right) ~\text{ for }~(\delta ,Q)\in (0,\infty )\\&\quad \times \left( 0,\displaystyle \frac{M_b}{\beta }\right) . \end{aligned}$$

A direct calculation gives $\partial u/\partial Q<0$. Note that $\displaystyle \lim _{Q\rightarrow 0^+}u(\delta ,Q)=(1+\gamma _1/M_b)(\delta ^{*}_0-\delta )>0$ if $0<\delta <\delta ^{*}_0$ and $\displaystyle \lim _{Q\rightarrow (M_b/\beta )^{-}}u(\delta ,Q)=-\infty$. Then for any fixed $\delta \in (0,\delta _0^*)$, there exists a unique positive $Q_\delta$ such that $u(\delta ,Q_\delta )=0$, and $A_2^*=Q_\delta /\delta$ is desired unique solution. This proves part (i).

Next we establish the stability of $E_2$. The stability of $E_2$ is determined by the eigenvalue problem

$$\begin{aligned}&\lambda \xi =-h_1\xi +h_2\zeta , \end{aligned}$$

(A.6a)

$$\begin{aligned}&\lambda \zeta = (-(1-p)\theta \delta +\theta h_1)\xi -\left( \displaystyle \frac{b}{x_e}+\theta h_2\right) \zeta +\displaystyle \frac{b}{x_e}\psi (0),\end{aligned}$$

(A.6b)

$$\begin{aligned}&\lambda \varphi =D_b\varphi ''(x)-v\varphi '(x)+ \left( h_3(x) -\delta \right) \varphi ,~0<x<x_h,\end{aligned}$$

(A.6c)

$$\begin{aligned}&\lambda \psi =D_m\psi ''(x)+ \left( p\theta \delta -\theta h_3(x)\right) \varphi ,~0<x<x_h,\end{aligned}$$

(A.6d)

$$\begin{aligned}&D_b\varphi '(0)-v\varphi (0)=D_b\varphi '(x_h)-v\varphi (x_h)=0,\end{aligned}$$

(A.6e)

$$\begin{aligned}&D_m\psi '(0)=b(\psi (0)-\zeta ),~\psi (x_h)=0, \end{aligned}$$

(A.6f)

where

$$\begin{aligned} \begin{aligned}&h_1=\displaystyle \frac{r\gamma _2\mu (A^*_2)N_2^*A^*_2}{N_2^*+\gamma _1},~ h_2=\displaystyle \frac{r\gamma _1A^*_2}{(N_2^*+\gamma _1)^2}\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_2^*))\mathrm{d}x,\\&h_3(x)=rf(M_2^*)g(I(x,A^*_2,0)) \end{aligned} \end{aligned}$$

(A.7)

and

$$\begin{aligned} \mu (A^*_2)=\displaystyle \frac{l\left( \ln \displaystyle \frac{I_{in}+\gamma _2}{I_{in}e^{-x_e(K_{bg}+lA_2^*)}+\gamma _2}-\displaystyle \frac{I_{in}x_e(K_{bg}+lA^*_2)e^{-x_e(K_{bg}+lA^*_2)}}{I_{in}e^{-x_e(K_{bg}+lA^*_2)}+\gamma _2} \right) }{x_e(K_{bg}+lA^*_2)^2}>0.\nonumber \\ \end{aligned}$$

(A.8)

Again let $\lambda _1$ be the eigenvalue of (A.6) with largest real part, and let $(\xi ,\zeta ,\varphi ,\psi )$ be the corresponding eigenfunction. Note that the eigenvalue problem (A.6) is partially decoupled. We consider the following two cases: (i) $\varphi \not \equiv 0$; (ii) $\varphi \equiv 0$.

Case (i): $\varphi \not \equiv 0$. Then, the stability of $E_2$ is determined by the subsystem

$$\begin{aligned} \begin{aligned}&\lambda \varphi =D_b\varphi ''(x)-v\varphi '(x)+\left( h_3(x)-\delta \right) \varphi ,~0<x<x_h,\\&D_b\varphi '(0)-v\varphi (0)=D_b\varphi '(x_h)-v\varphi (x_h)=0, \end{aligned} \end{aligned}$$

(A.9)

with $\xi ,\zeta ,\psi$ solved from (A.6a), (A.6b), (A.6d) and (A.6f). Then, $\lambda _1=\lambda _1(D_b,v,h_3(x)-\delta )=\lambda _1(D_b,v,h_3(x))-\delta$.

Case (ii): $\varphi \equiv 0$. Then, $(\xi ,\zeta ,\psi )$ satisfies

$$\begin{aligned} \begin{aligned}&\lambda \xi =-h_1\xi +h_2\zeta ,\\&\lambda \zeta = (-(1-p)\theta \delta +\theta h_1)\xi -\left( \displaystyle \frac{b}{x_e}+\theta h_2\right) \zeta +\displaystyle \frac{b}{x_e}\psi (0),\\&\lambda \psi =D_m\psi ''(x),~0<x<x_h,\\&D_m\psi '(0)=b(\psi (0)-\zeta ),~\psi (x_h)=0. \end{aligned} \end{aligned}$$

(A.10)

By using similar arguments as those in Theorem 3.1 and Lemma A.1, we conclude that all eigenvalues of (A.10) have negative real parts. Based on the analysis above, $E_2$ is locally asymptotically stable with respect to (2.1) if (3.11) holds, while $E_2$ is unstable if (3.12) holds.

It can be shown that $A_2^*$ is strictly decreasing in $\delta$ and $\displaystyle \lim _{\delta \rightarrow \delta _0^{*-}}A_2^*=0$. This also implies that $\displaystyle \lim _{\delta \rightarrow \delta _0^{*-}}N_2^*=M_b$ and $\displaystyle \lim _{\delta \rightarrow \delta _0^{*-}}M_2^*(x)=M_b$. Hence, when $\delta \rightarrow \delta _0^{*-}$, we have

$$\begin{aligned} \begin{aligned}&\lim _{\delta \rightarrow \delta _0^{*-}}\left[ \lambda _1\left( D_b,v,rf(M_2^*)g(I(x,A_2^*,0) \right) -\delta \right] \\ =&\lim _{\delta \rightarrow \delta _0^{*-}}\left[ \lambda _1\left( D_b,v,rf(M_2^*)g(I(x,A_2^*,0) \right) \right. \\&\left. -\lambda _1\left( D_b,v,f(N_2^*)\displaystyle \frac{r}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_2^*))\mathrm{d}x\right) \right] \\ =&\lambda _1\left( D_b,v,rf(M_b)g(I(x,0,0)\right) -\lambda _1\left( D_b,v,f(M_b)\displaystyle \frac{r}{x_e} \displaystyle \int _{-x_e}^0g(I(x,0))\mathrm{d}x\right) <0. \end{aligned} \end{aligned}$$

Therefore, (3.11) holds and $E_2$ is locally asymptotically stable when $\delta$ is sufficiently close to $\delta _0^*$. When $p=1$, $E_2\equiv (A_2^*,M_b,0,M_b)$ where $A_2^*$ satisfies (3.10) with $N_2^*=M_b$. From (3.5) and the monotonicity of the principal eigenvalue on the weight functions, we have

$$\begin{aligned} \lambda _1\left( D_b,v,rf(M_b)g(I(x,A_2^*,0)\right) < \lambda _1\left( D_b,v,f(M_b)\displaystyle \frac{r}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_2^*))\mathrm{d}x\right) =\delta , \end{aligned}$$

so (3.11) holds and $E_2$ is locally asymptotically stable for any $0<\delta <\delta _0^*$. $\square$

To obtain the global existence of $E_3$, we first establish the following a priori estimates for nonnegative solutions $(N_3,B_3(x),M_3(x))$ of (3.3).

Lemma A.2

Assume that $(N_3,B_3(x),M_3(x))\in {\mathbb {R}}^+\times Y^2$ is a nonnegative solution of (3.3) with $N_3>0,B_3,M_3\not \equiv 0$. Then,

• (i) $0<\delta <\delta _{**}$, where $\delta _{**}$ is defined in (3.5);

• (ii) For any $\epsilon >0$, there exists a positive constant $K(\epsilon )$ such that $0<B_3(x)\le K(\epsilon )$ and $0<N_3,M_3(x)\le M_b+\theta x_h^2(r+p\delta _{**})K(\epsilon )/D_m$ on $[0,x_h]$ for $\delta \in [\epsilon ,\delta _{**})$.

Proof

(i) Let $R(x)=B_3(x)e^{-(v/D_b)x}$. Then, R(x) satisfies

$$\begin{aligned} {\left\{ \begin{array}{ll} -D_bR''(x)-vR'(x)+\delta R=rf(M_3)g(I(x,0,B_3))R\ge 0,&{}0<x<x_h,\\ R'(0)=0,~R'(x_h)=0.&{} \end{array}\right. } \end{aligned}$$

It follows from the strong maximum principle that $R(x)>0$ on $[0,x_h]$, and consequently, $B_3(x)>0$ on $[0,x_h]$. From (3.3), we have

$$\begin{aligned} \lambda _1\left( D_b,v,rf(M_3)g(I(x,0,B_3))\right) =\delta \end{aligned}$$

with corresponding principal eigenfunction $B_3$. It follows from the monotonicity of principal eigenvalue with respect to the weight functions that

$$\begin{aligned} \delta =\lambda _1\left( D_b,v,rf(M_3)g(I(x,0,B_3))\right) <\lambda _1\left( D_b,v,rg(I(0,0,0))\right) =\delta _{**}. \end{aligned}$$

(ii) Fix $\epsilon >0$. If $B_3$ is not bounded for all $\delta \in [\epsilon ,\delta _{**})$, then there are a sequence $\delta ^i\in [\epsilon ,\delta _{**})$ and corresponding positive solutions $(N_3^i,B_3^i(x),M^i_3(x))$ such that $\Vert B_3^i\Vert _\infty \rightarrow \infty$ and $\delta ^i\rightarrow \delta ^0\in [\epsilon ,\delta _{**}]$ as $i\rightarrow \infty$. Let $\mu _i=B_3^i/\Vert B_3^i\Vert _\infty$. From the equation of $B_3$ in (3.3), we get

$$\begin{aligned} {\left\{ \begin{array}{ll} -D_b\mu _i''(x)+v\mu _i'(x)=[rf(M_3^i)g(I(x,0,B_3^i))-\delta ^i] \mu _i(x)=0,&{}0<x<x_h,\\ D_b\mu _i'(0)-v\mu _i(0)=D_b\mu _i'(x_h)-v\mu _i(x_h)=0.&{} \end{array}\right. }\nonumber \\ \end{aligned}$$

(A.11)

Since the right hand side of (A.11) is uniformly bounded, by using $L^p$ theory for elliptic operators and by passing to a subsequence, we obtain that $\mu _i\rightarrow \mu$ in $W^{2,p}([0,x_h)]$ (and also in $C^{1,\alpha }([0,x_h])$ from Sobolev’s embedding) as $i\rightarrow \infty$. Since $rf(M_3^i)g(I(x,0,B_3^i))$ is bounded in $L^{\infty }([0,x_h])$, we may assume that $rf(M_3^i)g(I(x,0,B_3^i))\rightarrow l^0$ weakly in $L^2([0,x_h])$. Hence, $\mu$ satisfies (in the weak sense)

$$\begin{aligned} {\left\{ \begin{array}{ll} D_b\mu ''(x)-v\mu '(x)+(l^0-\delta ^0) \mu (x)=0,&{}0<x<x_h,\\ D_b\mu '(0)-v\mu (0)=D_b\mu '(x_h)-v\mu (x_h)=0.&{} \end{array}\right. } \end{aligned}$$

(A.12)

From the strong maximum principle, $\mu (x)>0$ on $[0,x_h]$ since $\mu \ge 0$ and $\Vert \mu \Vert _\infty =1$. This implies that $B_3^i=||B_3^i||_{\infty }\mu _i\rightarrow \infty$ uniformly on $[0,x_h]$, and thus, $l^0=0$. Integrating (A.12) we obtain $0=\delta ^0\displaystyle \int _0^{x_h}\mu (x)\mathrm{d}x>0.$ That is a contradiction. Therefore, there is a positive constant $K(\epsilon )$ such that $0<B_3(x)\le K(\epsilon )$ on $[0,x_h]$ for all $\delta \in [\epsilon ,\delta _{**})$.

From the strong maximum principle, we have $M_3(x)>0$ on $[0,x_h]$. For any $x\in [0,x_h]$, we have

$$\begin{aligned} |D_mM'_3(x)|&=\left| D_m\displaystyle \int _0^{x}M''_3(s)ds\right| = \theta \left| \displaystyle \int _0^x\left( rf(M_3)g(I(x,0,B_3))-p\delta \right) B_3ds\right| \\&\le \theta x_h(r+p\delta _{**})K(\epsilon ) \end{aligned}$$

and consequently

$$\begin{aligned} \begin{aligned} |M_3(x)|&=|M_3(x_h)+M_3(x)-M_3(x_h)|\le |M_3(x_h)|+|M_3(x)-M_3(x_h)|\\&\le M_b+\left| \displaystyle \int _x^{x_h}M_3'(s)ds \right| \le M_b+\theta x_h^2(r+p\delta _{**})K(\epsilon )/D_m. \end{aligned} \end{aligned}$$

(A.13)

From $M_3(0)=N_3$, we have $0<N_3\le M_b+\theta x_h^2(r+p\delta _{**})K(\epsilon )/D_m$. $\square$

Proof of Theorem 3.5

We prove the existence of $E_3$ using bifurcation theory and show that the solution $E_3$ bifurcates from the line of nutrient-only semi-trivial steady state $E_1$ at $\delta =\delta _{*}$ with parameter $\delta$. We first use local bifurcation theory in Crandall and Rabinowitz (1971) to show that $E_3$ bifurcates from the line of $E_1$ at $\delta =\delta _{*}$.

Recall function spaces $X_1,X_2$ defined in (3.13) and define $Y=C([0,x_h])$. We define a nonlinear mapping $F:{\mathbb {R}}^+\times {\mathbb {R}}\times X_1\times X_2\rightarrow {\mathbb {R}}\times Y^2\times {\mathbb {R}}^2$ by

$$\begin{aligned} F(\delta ,N,B,M)=\begin{pmatrix} M(0)-N\\ D_bB''-vB'+rBf(M)g(I(x,0,B))-\delta B\\ D_mM''+\theta p\delta B-\theta rBf(M)g(I(x,0,B))\\ D_mM'(0)-b(M(0)-N)\\ M(x_h)-M_b \end{pmatrix}. \end{aligned}$$

It is easy to see that $F(\delta ,M_b,0,M_b)=0$. Let $H:=F_{(N,B,M)}(\delta _*,M_b,0,M_b)$ be the Frechét derivative of F with respect to (N, B, M) at $(\delta _*,M_b,0,M_b)$. For any $(\zeta ,\varphi ,\psi )\in {\mathbb {R}}\times X_1\times X_2$, we have

$$\begin{aligned} \begin{aligned} H[\zeta ,\varphi ,\psi ] =\begin{pmatrix}\psi (0)-\zeta \\ D_b\varphi ''(x)-v\varphi '(x)+\left( rf(M_b)g(I(x,0,0))-\delta _*\right) \varphi \\ D_m\psi ''(x)+\left( \theta p\delta _*-\theta rf(M_b)g(I(x,0,0))\right) \varphi \\ D_m\psi '(0)-b(\psi (0)-\zeta )\\ \psi (x_h) \end{pmatrix}. \end{aligned} \end{aligned}$$

(A.14)

If $(\zeta _1,\varphi _1,\psi _1)\in \ker H$, then

$$\begin{aligned}&\psi _1(0)-\zeta _1=0, \end{aligned}$$

(A.15a)

$$\begin{aligned}&D_b\varphi ''_1(x)-v\varphi '_1(x)+\left( rf(M_b)g(I(x,0,0))-\delta _*\right) \varphi _1=0,~0<x<x_h, \end{aligned}$$

(A.15b)

$$\begin{aligned}&D_m\psi ''_1(x)+\left( \theta p\delta _*-\theta rf(M_b)g(I(x,0,0))\right) \varphi _1=0,~0<x<x_h, \end{aligned}$$

(A.15c)

$$\begin{aligned}&\psi '_1(0)=0,~\psi _1(x_h)=0. \end{aligned}$$

(A.15d)

Recall that $\delta _*$ is the principal eigenvalue of (3.4) with $q(x)=rf(M_b)g(I(x,0,0))$, and let ${\bar{\varphi }}\in X_1$ be the corresponding positive eigenfunction for the principal eigenvalue $\delta _*$. Then, ${\bar{\varphi }}$ is the unique solution of (A.15b) up to a constant multiple. And there exist unique functions ${\bar{\psi }}\in X_2$ and ${\bar{\zeta }}\in {\mathbb {R}}$ satisfying (A.15c), (A.15d) and (A.15a). Hence, $\dim \ker H=1$ and $\ker H={{\,\mathrm{span}\,}}\{({\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }})\}.$

If $(\sigma _1,\sigma _2,\sigma _3,\sigma _4,\sigma _5)\in {{\,\mathrm{range}\,}}H$, then there exists $(\zeta _2,\varphi _2,\psi _2)\in {\mathbb {R}}\times X_1\times X_2$ such that

$$\begin{aligned} \begin{aligned}&\psi _2(0)-\zeta _2=\sigma _1,\\&D_b\varphi ''_2(x)-v\varphi '_2(x)+\left( rf(M_b)g(I(x,0,0))-\delta _*\right) \varphi _2=\sigma _2(x),~0<x<x_h,\\&D_m\psi ''_2(x)+\left( \theta p\delta _*-\theta rf(M_b)g(I(x,0,0))\right) \varphi _2=\sigma _3(x),~0<x<x_h,\\&D_m\psi '_2(0)-b(\psi _2(0)-\zeta _2)=\sigma _4,~\psi _2(x_h)=\sigma _5. \end{aligned} \end{aligned}$$

(A.16)

Multiplying both sides of (A.15b) and the second equation of (A.16) by $\varphi _2e^{-(v/D_b)x}$ and $\varphi _1e^{-(v/D_b)x}$, respectively, subtracting, and integrating on $[0,x_h]$, we obtain

$$\begin{aligned}&\displaystyle \int _0^{x_h}\sigma _2(x)e^{-(v/D_b)x}\varphi _1(x)\mathrm{d}x\\=&D_b\displaystyle \int _0^{x_h}\left( \left( \varphi '_2(x) e^{-(v/D_b)x}\right) '\varphi _1(x)-\left( \varphi '_1(x) e^{-(v/D_b)x}\right) '\varphi _2(x)\right) \mathrm{d}x\\ =&D_b\varphi '_2(x)e^{-(v/D_b)x}\varphi _1(x)\Big |_0^{x_h}-D_b\varphi '_1(x)e^{-(v/D_b)x} \varphi _2(x)\Big |_0^{x_h}=0. \end{aligned}$$

This implies that

$$\begin{aligned} \displaystyle \int _0^{x_h}\sigma _2(x)e^{-(v/D_b)x}{\bar{\varphi }}(x)\mathrm{d}x=0 \end{aligned}$$

and

$$\begin{aligned} {{\,\mathrm{range}\,}}H=\left\{ (\sigma _1,\sigma _2,\sigma _3,\sigma _4,\sigma _5)\in {\mathbb {R}}\times Y^2\times {\mathbb {R}}^2:\displaystyle \int _0^{x_h}\sigma _2(x)e^{-(v/D_b)x}{\bar{\varphi }}(x)\mathrm{d}x=0\right\} . \end{aligned}$$

Then, ${{\,\mathrm{codim}\,}}{{\,\mathrm{range}\,}}H=1$. Moreover, we also have

$$\begin{aligned} F_{\delta (N,B,M)}(\delta _*,M_b,0,M_b) ({\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }})=(0,-{\bar{\varphi }},p\theta {\bar{\varphi }},0,0), \end{aligned}$$

so $F_{\delta ,(N,B,M)}(\delta _*,M_b,0,M_b)({\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }}) \not \in {{\,\mathrm{range}\,}}H$ as $\displaystyle \displaystyle \int _0^{x_h}\sigma _2(x)e^{-(v/D_b)x}{\bar{\varphi }}^2(x)\mathrm{d}x\ne 0$.

From Theorem 1.7 in Crandall and Rabinowitz (1971), there exists a positive constant $\varepsilon _3>0$ such that all solutions of (3.3) near $(\delta _*,M_b,0,M_b)$ lie on a smooth curve

$$\begin{aligned} \Gamma _3=\{(\delta _3(s),N_3^*(s),B_3^*(s,x),M_3^*(s,x)):0<s<\varepsilon _3\} \end{aligned}$$

with the form

$$\begin{aligned} {\left\{ \begin{array}{ll} &{}N_3^*(s)=M_b+s{\bar{\zeta }}+o(s),\;\; B_3^*(s,x)=s{\bar{\varphi }}(x)+o(s),\\ &{}M_3^*(s,x)=M_b+s{\bar{\psi }}(x)+o(s). \end{array}\right. } \end{aligned}$$

(A.17)

Let ${\bar{l}}$ be a linear functional on ${\mathbb {R}}\times Y^2\times {\mathbb {R}}^2$ by

$$\begin{aligned} \langle {\bar{l}},(\sigma _1,\sigma _2,\sigma _3,\sigma _4,\sigma _5)\rangle = \displaystyle \int _0^{x_h}\sigma _2(x)e^{-(v/D_b)x}{\bar{\varphi }}(x)\mathrm{d}x. \end{aligned}$$

Then, from Liu et al. (2007), we have

$$\begin{aligned} \begin{aligned} \delta _3'(0)&=-\displaystyle \frac{\left\langle {\bar{l}},F_{(N,B,M)(N,B,M)}\left( \delta _*, M_b,0,M_b\right) [{\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }}]^2\right\rangle }{2\left\langle {\bar{l}}, F_{\delta ,(N,B,M)}\left( \delta _*, M_b,0,M_b\right) [{\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }}]\right\rangle }\\&=-\displaystyle \frac{\displaystyle \int _0^{x_h}\displaystyle \frac{rI(x,0,0)e^{-(v/D_b)x}}{(M_b+\gamma _1)(I(x,0,0)+\gamma _2)} \mu (x,{\bar{\varphi }},{\bar{\psi }})\mathrm{d}x}{\displaystyle \int _0^{x_h}e^{-(v/D_b)x} {\bar{\varphi }}^2(x)\mathrm{d}x}, \end{aligned} \end{aligned}$$

(A.18)

where

$$\begin{aligned} \mu (x,{\bar{\varphi }},{\bar{\psi }})= \displaystyle \frac{\gamma _2 lM_b}{I(x,0,0)+\gamma _2}{\bar{\varphi }}^2(x)\displaystyle \int _0^x {\bar{\varphi }}(s)ds- \displaystyle \frac{\gamma _1}{M_b+\gamma _1}{\bar{\varphi }}^2(x){\bar{\psi }}(x),~x\in [0,x_h]. \end{aligned}$$

We claim that ${\bar{\psi }}(x)\le 0$ on $x\in [0,x_h]$. Let $y(x)=rf(M_b)g(I(x,0,0)).$ Then y(x) is a strictly decreasing function on $x\in [0,x_h]$. It follows from (A.15b)-(A.15d) that ${\bar{\psi }}'(x_h)\ge 0$. From (A.15c), we have ${\bar{\psi }}''(x)>0$ on $x\in [0,x_h]$ if $p\le y(x_h)/\delta _*$. Combining with its boundary conditions ${\bar{\psi }}'(x)=0,{\bar{\psi }}(x_h)=0$, we get ${\bar{\psi }}(x)\le 0$ on $x\in [0,x_h]$. If $p>y(x_h)/\delta _*$, then there is a unique $x_0\in (0,x_h)$ such that ${\bar{\psi }}''(x)>0$ on $x\in (0,x_0)$ and ${\bar{\psi }}''(x)<0$ on $x\in (x_0,x_h)$. This means that ${\bar{\psi }}'(x)$ is a strictly increasing function on $x\in [0,x_0)$ and a strictly decreasing function on $x\in (x_0,x_h]$. Hence ${\bar{\psi }}'(x)>0$ on $x\in (0,x_h)$ since ${\bar{\psi }}'(0)=0$ and ${\bar{\psi }}'(x_h)\ge 0$. It follows from ${\bar{\psi }}(x_h)=0$ that ${\bar{\psi }}(x)\le 0$ on $x\in [0,x_h]$. From (A.18), we have $\delta _3'(0)<0$, which implies that the bifurcation at $(\delta _*,M_b,0,M_b)$ is backward. This completes the proof of part (iii).

Now we turn to global bifurcation of solutions of (3.3). By using Theorem 3.3 and Remark 3.4 in Shi and Wang (2009), we conclude that there exists a connected component $\Lambda ^+$ of $\Lambda$ containing $\Gamma _3$, and the closure of $\Lambda ^+$ includes the bifurcation point $(\delta _*,M_b,0,M_b)$. Moreover, $\Lambda ^+$ satisfies one of the following three alternatives:

• (a) it is not compact in ${\mathbb {R}}^+\times {\mathbb {R}}\times X_1\times X_2$;

• (b) it contains another bifurcation point $({\hat{\delta }},M_b,0,M_b)$ with ${\hat{\delta }}\ne \delta _*$;

• (c) it contains a point $(\delta ,M_b+N_3^{**},B_3^{**}(x), M_b+M_3^{**}(x))$ with $0\ne (N_3^{**},B_3^{**}(x),M_3^{**}(x))\in {\mathcal {W}}$, where ${\mathcal {W}}$ is a closed complement of $\ker H={{\,\mathrm{span}\,}}({\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }})$ in ${\mathbb {R}}\times X_1\times X_2$.

If the alternative (c) holds, then

$$\begin{aligned} \int _0^{x_h}B_3^{**}(x){\bar{\varphi }}(x)\mathrm{d}x=0. \end{aligned}$$

(A.19)

But it follows from Lemma A.2 that $B_3^{**}(x)>0$ on $[0,x_h]$, and ${\bar{\varphi }}(x)>0$ as it is a positive eigenfunction. This is a contradiction to (A.19), which means that (c) cannot happen. Suppose that the alternative (b) occurs, then (A.15) has a nonzero solution $({{\hat{\zeta }}},{{\hat{\varphi }}},{{\hat{\psi }}})$ with $\delta _*$ replaced by ${\hat{\delta }}$, ${\hat{\varphi \in }} X_1$, and ${{\hat{\varphi }}}>0$. But the eigenvalue problem (3.4) has only one eigenvalue with positive eigenfunction, hence (b) cannot happen either.

Therefore, the alternative (a) must happen, and $\Lambda ^+$ is not compact in ${\mathbb {R}}^+\times {\mathbb {R}}\times X_1\times X_2$. From Lemma A.2, $(N_3^*,B_3^*(x),M_3^*(x))$ is bounded on $[0,x_h]$ for $\delta \in [\epsilon ,\delta _{**})$ with any $\epsilon >0$, and (3.3) has no nonnegative solution when $\delta >\delta _{**}$. Hence the projection of $\Lambda ^+$ onto $\delta$-axis is contained in $(0,\delta _{**})$, but contains $[\epsilon ,\delta _*)$ for any $\epsilon >0$, so the projection of $\Lambda ^+$ onto $\delta$-axis contains $(0,\delta _*)$. This also implies that there is at least one positive solution of (3.3) on $\Lambda ^+$ for any $\delta \in (0,\delta _*)$. This completes the proof of parts (i) and (ii).

$\square$

To prove the existence of coexistence steady state $E_4$, we first prove the following elementary result.

Lemma A.3

Suppose that $0<a<x_e(\delta _0^*-\delta _*)$ holds, then

$$\begin{aligned} {\left\{ \begin{array}{ll} D_b\varphi ''(x)-v\varphi '(x)+\left( rf(M_b)g(I(x,0,0))-\delta _a^*\right) \varphi =0,&{}~0<x<x_h,\\ D_b\varphi '(0)-v\varphi (0)=-a,~D_b\varphi '(x_h)-v\varphi (x_h)=0.&{} \end{array}\right. }\nonumber \\ \end{aligned}$$

(A.20)

has a unique positive solution ${\hat{\varphi }}(x)$.

Proof

By using the transform $\phi (x)=\varphi (x)e^{-(v/D_b)x}$, we get

$$\begin{aligned} {\left\{ \begin{array}{ll} D_b\phi ''(x)+v\phi '(x)+\left( rf(M_b)g(I(x,0,0))-\delta _a^*\right) \phi =0,&{}~0<x<x_h,\\ -D_b\phi '(0)-a=0,~\phi '(x_h)=0.&{} \end{array}\right. }\nonumber \\ \end{aligned}$$

(A.21)

It is clear that 0 is a lower solution of (A.21). Let ${\tilde{\varphi }}(x)={\tilde{\phi }}(x)e^{(v/D_b)x}$ be the positive eigenfunction of (3.4) corresponding to $\lambda _1(D_b,v,rf(M_b)g(I(x,0,0))=\delta _*$ with $\max _{x\in [0,x_h]}{\tilde{\phi }}(x)=1$. Define ${\bar{\phi }}(x)=K [{\tilde{\phi }}(x)+\epsilon (x-x_h)^2]$, where K and $\epsilon$ are positive constants to be specified. Then

$$\begin{aligned} \begin{aligned}&D_b{\bar{\phi }}''(x)+v{\bar{\phi }}'(x)+\left( rf(M_b)g(I(x,0,0)) -\delta _a^*\right) {\bar{\phi }}\\ =&(\delta _*-\delta _a^*)K{\tilde{\phi }}+\epsilon K \left[ 2D_b+2v(x-x_h)+\left( rf(M_b)g(I(x,0,0)) -\delta _a^*\right) (x-x_h)^2\right] \\\le&0, \;\;\;\; 0<x<x_h, \end{aligned} \end{aligned}$$

if $\epsilon >0$ is chosen sufficiently small since $\delta _*-\delta _a^*<0$ and ${\tilde{\phi }}>0$. We can further choose $K>0$ so that $-D_b{\bar{\phi }}'(0)-a=2D_bK\epsilon x_h-a=0$ and ${\bar{\phi }}'(x_h)=0.$

Hence ${\bar{\phi }}$ is an upper solution of (A.21). From Theorem 3.2.1 in Pao (1992), there is a solution ${\hat{\phi }}$ of (A.21) satisfying $0\le {\hat{\phi }}\le {\bar{\phi }}$. It follows from the strong maximum principle that ${\hat{\phi }}>0$. Note that (A.21) is a linear ODE with nonhomogeneous boundary conditions. Hence ${\hat{\phi }}$ is unique. This implies that there exists a unique positive solution ${\hat{\varphi }}={\hat{\phi }}e^{(v/D_b)x}$ of (A.20). $\square$

To obtain the global existence of $E_4$, we first establish the following a priori estimates for nonnegative solutions $(A_4,N_4,B_4(x),M_4(x))$ of (3.14).

Lemma A.4

Assume that $(A_4,N_4,B_4,M_4)\in ({\mathbb {R}}^+)^2\times Y^2$ is a nonnegative solution of (3.14) with $A_4,N_4>0$ and $B_4,M_4\not \equiv 0$. Then

• (i) $0<\delta <\delta _a^{**}$, where $\delta _a^{**}$ is defined in (3.5);

• (ii) $0<A_4<K$, where K satisfies $\displaystyle \frac{r}{x_e}\int _{-x_e}^0g(I(x,K))\mathrm{d}x=\delta +\displaystyle \frac{a}{x_e}$;

• (iii) For any $\epsilon >0$, there exists a positive constant $C(\epsilon )>0$ such that

  $$\begin{aligned}&0<B_4(x)\le C(\epsilon ),\\&0<N_4,M_4(x)<M_b+x_h^2\theta (r+p\delta _a^{**})C(\epsilon )/D_m+ x_h\theta (x_e(1-p)\delta _a^{**}+a)K/D_m \end{aligned}$$

  on $[0,x_h]$ for $\delta \in [\epsilon ,\delta _a^{**})$.

Proof

It follows from the first equation of (3.14) that (i) and (ii) hold. By applying similar arguments to those in Lemma A.2, we conclude that $M_4>0$ and for any $\epsilon >0$, there exists a positive constant $C(\epsilon )>0$ such that $0<B_4\le C(\epsilon )$. From the second and fourth equations of (3.14), we have

$$\begin{aligned} b(M_4(0)-N_4)=\theta ( x_e(1-p)\delta +a)A_4 \end{aligned}$$

and

$$\begin{aligned} |D_mM'_4(x)|&=\left| D_m\displaystyle \int _0^{x}M''_4(s)ds+D_mM_4'(0)\right| \\&\le \theta \left| \displaystyle \int _0^x\left( rf(M_4)g(I(x,A_4,B_4))-p\delta \right) B_4ds\right| + b|M_4(0)-N_4|\\&\le x_h\theta (r+p\delta _a^{**})C(\epsilon )+\theta (x_e(1-p)\delta _a^{**}+a)K. \end{aligned}$$

By (A.13), we conclude that (iii) holds. $\square$

Proof of Theorem 3.7

Define a nonlinear mapping $G:{\mathbb {R}}^+\times {\mathbb {R}}^2\times X_3\times X_2\rightarrow {\mathbb {R}}^2\times Y^2\times {\mathbb {R}}^3$ by

$$\begin{aligned} G(\delta ,A,N,B,M)=\begin{pmatrix} \left( rf(N)\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A))\mathrm{d}x-\delta -\displaystyle \frac{a}{x_e}\right) A\\ \displaystyle \frac{b}{x_e}(M(0)-N)+\theta \left( p \delta -rf(N)\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A))\mathrm{d}x\right) A\\ D_bB''-vB'+rBf(M)g(I(x,A,B))-\delta B\\ D_mM''+\theta \left( p\delta - rf(M)g(I(x,A,B))\right) B\\ D_bB'(0)-vB(0)+aA\\ D_mM'(0)-b(M(0)-N)\\ M(x_h)-M_b \end{pmatrix},\nonumber \\ \end{aligned}$$

(A.22)

where $X_3:=\{u\in C^2([0,x_h]):D_bu'(x_h)-vu(x_h)=0\}$, and $X_2,Y$ are defined in (3.13). It follows that $G(\delta ,0,M_b,0,M_b)=0$. Let $L:=G_{(A,N,B,M)}(\delta _a^*,0,M_b,0,M_b)$ be the Frechét derivative of G with respect to (A, N, B, M). For any $(\xi ,\zeta ,\varphi ,\psi )\in {\mathbb {R}}^2\times X_3\times X_2$, we have

$$\begin{aligned} \begin{aligned} L[\xi ,\zeta ,\varphi ,\psi ] =\begin{pmatrix}0\\ -\theta \left( (1-p)\delta _a^*+\displaystyle \frac{a}{x_e}\right) \xi -\displaystyle \frac{b}{x_e}\zeta +\displaystyle \frac{b}{x_e}\psi (0)\\ D_b\varphi ''(x)-v\varphi '(x)+(rf(M_b)g(I(x,0,0))-\delta _a^*)\varphi \\ D_m\psi ''(x)+\theta (p\delta _a^*-rf(M_b)g(I(x,0,0)))\varphi \\ D_b\varphi '(0)-v\varphi (0)+a\xi \\ D_m\psi '(0)-b(\psi (0)-\zeta )\\ \psi (x_h) \end{pmatrix}. \end{aligned} \end{aligned}$$

If $(\xi _1,\zeta _1,\varphi _1,\psi _1)\in \ker L$, then

$$\begin{aligned} \begin{aligned}&-\theta \left( (1-p)\delta _a^*+\displaystyle \frac{a}{x_e}\right) \xi _1-\displaystyle \frac{b}{x_e}\zeta _1 +\displaystyle \frac{b}{x_e}\psi _1(0)=0,\\&D_b\varphi ''_1(x)-v\varphi '_1(x)+(rf(M_b)g(I(x,0,0))-\delta _a^*)\varphi _1=0,~0<x<x_h,\\&D_m\psi ''_1(x)+\theta (p\delta _a^*-rf(M_b)g(I(x,0,0)))\varphi _1=0,~0<x<x_h,\\&D_b\varphi '_1(0)-v\varphi _1(0)+a\xi _1=0,~D_m\psi '_1(0)-b(\psi _1(0)-\zeta _1)=0,~\psi _1(x_h)=0. \end{aligned} \end{aligned}$$

(A.23)

Let $\xi _1=1$. Then, from Lemma A.3, $\varphi _1={\hat{\varphi }}>0$ can be uniquely solved, so are $\zeta _1={\hat{\zeta }}$ and $\psi _1={\hat{\psi }}$. Similar to those in the proof of (i) in Theorem 3.5, we have ${\hat{\psi }}\le 0$ and ${\hat{\zeta }}<0$. If $\xi _1=0$, we can conclude that $\varphi _1\equiv 0$ as $\delta _*<\delta _a^*$ and consequently $\zeta _1=\psi _1=0$. Hence, $\dim \ker L=1$ and $\ker L={{\,\mathrm{span}\,}}\{(1,{\hat{\zeta }},{\hat{\varphi }},{\hat{\psi }})\}$. It is also easy to observe that ${{\,\mathrm{codim}\,}}{{\,\mathrm{range}\,}}L=1$ as

$$\begin{aligned} {{\,\mathrm{range}\,}}L=&\left\{ (\rho _1,\rho _2,\rho _3,\rho _4,\rho _5,\rho _6,\rho _7)\in {\mathbb {R}}^2\times Y^2\times {\mathbb {R}}^3:\rho _1=0\right\} , \end{aligned}$$

and

$$\begin{aligned} G_{\delta (A,N,B,M)}(\delta _a^*,0,M_b,0,M_b) (1,{\hat{\zeta }},{\hat{\varphi }},{\hat{\psi }}) =(-1,\theta p,-{\hat{\varphi }},\theta p{\hat{\varphi }},0,0,0)\not \in {{\,\mathrm{range}\,}}L. \end{aligned}$$

By applying Theorem 1.7 in Crandall and Rabinowitz (1971), there exists a positive constant $\varepsilon _4>0$ such that all solutions of (3.14) near $(\delta _*,M_b,0,M_b)$ lie on a smooth curve

$$\begin{aligned} \Gamma _4=\{(\delta _4(s),A_4^*(s),N_4^*(s),B_4^*(s,x),M_4^*(s,x)):0<s<\varepsilon _4\} \end{aligned}$$

with the form

$$\begin{aligned} {\left\{ \begin{array}{ll} &{}A_4^*(s)=s+o(s^2),\;\; N_4^*(s)=M_b+s{\hat{\zeta }}+o(s^2),\\ &{}B_4^*(s,x)=s{\hat{\varphi }}(x)+o(s^2),\;\; M_4^*(s,x)=M_b+s{\hat{\psi }}(x)+o(s^2). \end{array}\right. } \end{aligned}$$

(A.24)

Again the direction of bifurcation can be calculated by

$$\begin{aligned} \begin{aligned} \delta _4'(0)&=-\displaystyle \frac{\left\langle {\hat{l}},G_{(A,N,B,M)(A,N,B,M)}(\delta _a^*,0,M_b,0,M_b)[1,{\hat{\zeta }}, {\hat{\varphi }},{\hat{\psi }}]^2\right\rangle }{2\left\langle {\hat{l}}, G_{\delta (A,N,B,M)}(\delta _a^*,0,M_b,0,M_b)[1,{\hat{\zeta }}, {\hat{\varphi }},{\hat{\psi }}]\right\rangle }\\&=-\displaystyle \frac{r\gamma _2\mu (0)M_b}{M_b+\gamma _1}+\displaystyle \frac{r\gamma _1{\hat{\zeta }}}{(M_b+\gamma _1)^2}\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,0))\mathrm{d}x<0 \end{aligned} \end{aligned}$$

where ${\hat{l}}$ is a linear functional on ${\mathbb {R}}^2\times Y^2\times {\mathbb {R}}^3$ defined as $\langle {\hat{l}},(\rho _1,\rho _2,\rho _3,\rho _4,\rho _5,\rho _6,\rho _7)\rangle = \rho _1$ and $\mu (0)$ can been found in (A.8). This shows that the bifurcation at $(\delta _a^*,0,M_b,0,M_b)$ is backward. This proves part (iii). The proof of part (i) and (ii) is similar to the one for Theorem 3.5, so we omit it here. $\square$

Appendix B

We briefly describe the numerical algorithm used in the paper. Divide the interval $[0,x_h]$ to n equal size subintervals $[x^i,x^{i+1}]$, $i=1,2,\ldots ,n,n+1$ with $x^i=x^{i-1}+\Delta x, x^1=0$ and $\Delta x=x_h/n$ (grid size). We also denote the step size (in the time direction) by $\Delta t$. Let

$$\begin{aligned} u_j^i=u(x^i,t_j),~t_j=t_{j-1}+\Delta t,~x^i=x^{i-1}+\Delta x,~1\le j\le m,~1\le i\le n+1. \end{aligned}$$

We also define $u^0_j$ to be the value at $x^0=-\Delta x$ and $u^{n+2}_j$ to be the value at $x^{n+2}=x_h+\Delta x$. We use backward differences and central differences as

$$\begin{aligned} u_t(x^i,t_j)=\displaystyle \frac{u_{j}^i-u_{j-1}^i}{\Delta t},~u_x(x^i,t_j)=\displaystyle \frac{u_{j}^i-u_{j}^{i-1}}{\Delta x},~ u_{xx}(x^i,t_j)=\displaystyle \frac{u_{j}^{i+1}-2u_{j}^i+u_j^{i-1}}{(\Delta x)^2}. \end{aligned}$$

Then, model (2.1) becomes

$$\begin{aligned} \begin{aligned}&\displaystyle \frac{A_j-A_{j-1}}{\Delta t}=rA_jf(N_j)\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_j))\mathrm{d}x -\displaystyle \frac{a}{x_e}A_j-\delta A_j,\\&\displaystyle \frac{N_j-N_{j-1}}{\Delta t}=\displaystyle \frac{b}{x_e}(M^1_j-N_j)+\theta p\delta A_j-\theta rA_jf(N_j)\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_j))\mathrm{d}x,\\&\displaystyle \frac{B_{j}^i-B_{j-1}^i}{\Delta t}=D_b\displaystyle \frac{B_{j}^{i+1}-2B_{j}^i+B_j^{i-1}}{(\Delta x)^2}-v\displaystyle \frac{B_{j}^i-B_{j}^{i-1}}{\Delta x}\\&\qquad \qquad \qquad \quad +rB_j^if(M_j^i)g(I(x^i,A_j,B_j^i))-\delta B_j^i,\\&\displaystyle \frac{M_{j}^i-M_{j-1}^i}{\Delta t}=D_m\displaystyle \frac{M_{j}^{i+1}-2M_{j}^i+M_j^{i-1}}{(\Delta x)^2} +\theta p\delta B_j^i-\theta rB_j^if(M_j^i)g(I(x^i,A_j,B_j^i)),\\&D_b\displaystyle \frac{B_{j}^1-B_{j}^{0}}{\Delta x}-vB_j^1=-aA_j,~D_b\displaystyle \frac{B_{j}^{n+2}-B_{j}^{n+1}}{\Delta x}-vB_j^{n+1}=0,\\&D_m\displaystyle \frac{M_{j}^1-M_{j}^{0}}{\Delta x}=b(M_{j}^1-N_j),~D_m\displaystyle \frac{M_{j}^{n+2}-M_{j}^{n+1}}{\Delta x}=M_b,\\&\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_j))\mathrm{d}x=\frac{1}{x_e(K_{bg}+lA_j)}\ln \frac{I_{in}+\gamma _2}{I_{in}\exp (-x_e(K_{bg}+lA_j))+\gamma _2}, \\&I(x^i,A_j,B_j^i)=I_{in}\exp \left( -K_{bg}(x^i+x_e)-x_elA_j-l\sum _{k=0}^i B_j^k\right) , \end{aligned} \end{aligned}$$

for $1\le j\le m,~1\le i\le n+1$. The simulations in the paper are based on the above numerical algorithm and are implemented in MATLAB.