STATIS multivariate three-way method for evaluating quality of life after corneal surgery: Methodology and case study in Costa Rica

1.   Introduction

Plant-herbivore interactions have been studied previously by a huge number of investigators including mathematicians and ecologists by using differential equation models and theories in dynamical systems. One of the most commonly use differential equations used the traditional Holling Type II functional response, which assumes that the growth rate of herbivore is a monotonically increasing function of plant density. However, this will not be reasonable if the chemical defense of plants is considered, in which case the negative effect of plant toxin on herbivore can lead to a decrease in the growth rate when the plant density is high. To explore the impact of plant toxicity on the dynamics of plant-herbivore interactions, models that include a toxin-determined functional response are proposed. The toxin-determined functional response is a modification of the traditional Holling Type II response by including the negative effect of toxin on herbivore growth, which can overwhelm the positive effect of biomass ingestion at sufficiently high plant toxicant concentrations.

In this paper, we consider the following system of plant-herbivore interactions with toxin-determined functional responses:

where $\Omega$ is a bounded open domain in ${\mathbf R}^N$ ($N\geq1$) with a smooth boundary $\partial\Omega$; $U$ and $V$ are the densities of the plant and herbivore for $(x, \tau)\in\Omega\times(0, \infty)$ respectively; $D_1$ and $D_2$ stand for the diffusion rates of the plant (for seed dispersal) and the predator respectively. $A$ and $K$ are the intrinsic growth rate and carrying capacity of the plant species respectively; $e$ is the rate of encounter per unit plant per herbivore, $h$ is the average time required for handling one unit of plant biomass; $C$ is the conversion rate of the consumed plant biomass into new herbivore biomass; $E$ is the per capita death rate of the herbivore due to causes unrelated to plant toxicity. The term $1-eBU/(1+heU)$, captures the negative toxin effects, where it is required that $h/4 < B < h$.

System (1.1) has been studied extensively by several authors, but most of the research focuses either on the corresponding ODEs system or on the PDEs system but concentrating on the traveling wave solutions. For example, Feng, Liu and DeAngelis ^[1] studied the stability of the equilibrium solutions and Hopf bifurcations of the ODEs system; Liu, Feng, Zhu and DeAngelis ^[2] further performed detailed bifurcation analysis of the ODEs system reveals a rich array of possible behaviors including cyclical dynamics through Hopf bifurcations and homoclinic bifurcation; Castillo-Chavez, Feng and Huang ^[3] studied the global dynamics of the corresponding ODEs system. They were able to find necessary and sufficient condition on the nonexistence of a closed orbit via the transformation of the model to a new equivalent system. The Poincare-Bendixson theorem was used to show that the existence of a unique interior equilibrium point guarantees its global asymptotical stability whenever it is locally asymptotically stable. Zhao, Feng, Zheng and Cen ^[4] studied the existence of limit cycles and homoclinic bifurcation in the ODEs. By using the theory of rotated vector fields and the extended planar termination principle, they showed the existence of limit cycles and homo-clinic loop. Li, Feng, Swihart, Bryant and Huntly ^[5] considered an ODE system of $n$ plant species and one herbivore population, which exhibits a rich variety of complex dynamics including Hopf bifurcation and period-doubling bifurcations. Feng, Huang and DeAngelis ^[6] studied the reaction-diffusion models for plant-herbivore interactions with toxin-determined functional response with two plant species that have different levels of toxicity. It was shown that under suitable conditions, the system might have traveling wave solutions; More recently, Xiang, Wu and Wan ^[7] use the steady state bifurcation theory and Hopf bifurcation theory to study the existence of multiple bifurcations.

In this paper, however, we shall limit our attention to the occurrence of steady state solutions bifurcating from the positive constant equilibrium solution by using the classical local and global steady state bifurcation theory in the sense of Crandall and Rabinowitz (^[8,9]). By using the diffusion rate $d_2$ of the predator as the main bifurcation parameter, we are able to show the existence of global bifurcating branches consisting of positive non-constant steady state solutions, that is., for any $d_2$ larger than certain critical value, then the diffusive system will have at least a positive non-constant steady state solution. On the other hand, compared with the diffusion rate $d_1$ of the plant, the diffusion rate $d_2$ of the herbivore tend to be sufficiently large from practical point of view. Thus, for any fixed $d_1$, the dynamics of the system for sufficient larger $d_2$ will be of great concern. Note that one way to solve the steady states of (1.2) for the larger $d_2$ case is to use the "shadow system"approach ^[10], which basically considers the limiting system when the diffusion rate $d_2$ tends to infinity. The shadow system is presumably easier to analyze, and the wish is that properties of solutions to the shadow system reflect that of the solutions to the original system, at least for $d_2$ sufficiently large. To a certain extent, this is true. From this point of view, the shadow system approach is the main advantage of our work compared with other existing literatures when we consider the case of $d_2$ sufficiently large.

Based on the aforementioned observations, we are to prove the existence of the positive non-constant steady state solutions of both original system and its limiting system-the shadow system. To make our discussions simpler, we introduce the following change of variables:

Then, we can reduce system (1.1) to the following dimensionless form

where $\gamma\in(1/4, 1)$. System (1.2) is similar to the classical predator-prey system with Holling type-II functional responses where $\gamma = 0$. (For more details, we refer to ^[11,12] and references therein). To facilitate the needs to perform our bifurcation analysis, without loss of generality, we treat only the simpler case when $\Omega: = (0, \ell)$ for $\ell > 0$.

The remaining parts of this paper are organized as follows. In Section $2$, we consider the spatial patterns of the original system by performing detailed local and global steady state bifurcation to the original system; In Section 3, we study the spatial patterns of the shadow system (the limiting system when $d_2\rightarrow \infty$) by using the time-mapping methods; In Section 4, we end up our discussions by drawing some conclusions.

2.   The spatial patterns of the original system: Steady state bifurcation analysis

In this section, we shall perform detailed global steady state bifurcation analysis to the system by using $d_2$ as the bifurcation parameter. Without loss of generality, we shall limit our attention to the case of $\Omega = (0, \ell)$ for $\ell > 0$. The steady state problem of system (1.1) is governed by the following elliptic equations:

One can check that: if $1/4 < \gamma\leq1/2$ holds, then $\dfrac{u}{1+u}(1- \frac{\gamma u}{1+u})$ is increasing in $u$, while it is unimodal if $1/2 < \gamma < 1$ holds, more precisely, if $0 < u < 1/(2\gamma-1)$, it is increasing, while if $u > 1/(2\gamma-1)$, it is decreasing.

Clearly, if $\frac{1}{4} < \gamma\leq\dfrac{1}{2}$ and $\frac{\theta}{m}\in\bigg(0, \frac{k}{1+k}\bigg(1- \frac{\gamma k}{1+k}\bigg)\bigg)$ holds, then system (2.1) has a unique positive constant solution, denoted by $(\lambda, v_ \lambda)$, satisfying $0 < \lambda < k$ and

Define $X = \{(u, v): u, v\in C^2[0, \ell], u' = v' = 0 \text{ on}\; x = 0, \ell\}$ and $Y = L^2(0, \ell)\times L^2(0, \ell)$ be the Hilbert space with the inner product

for $U_1: = (u_1, v_1)\in Y, U_2: = (u_2, v_2)\in Y$. Define the mapping $\mathcal F: (0, \infty)\times X\rightarrow Y$ by

Clearly, for all $d_2 > 0$, we have $\mathcal F(d_2, \lambda, v_ \lambda) = 0$. The linearized operator evaluated at $(\lambda, v_ \lambda)$ is given by

where $\Delta: = \frac{\partial^2}{\partial x^2}$ and

By ^[7,10,12], the eigenvalues of $L$ can be determined by the eigenvalues of $L_n$ defined by

where $\tau_n$ is the eigenvalue of $-\Delta$ on $\Omega$ subject to the homogenous Neumann boundary conditions satisfying $0 = \tau_0 < \tau_1 < \tau_2 < \cdots$. The characteristic equation of $L_n$ is given by

where

One can check that: if $k\leq1/(1+\gamma)$ holds, then for all $\lambda\in(0, k)$, $\sigma_1(\lambda) < 0$; while if $k > 1/(1+\mu)$ holds, then there exists a $\widehat \lambda\in(0, k)$, such that for any $\lambda\in(0, \widehat \lambda)$, $\sigma_1(\lambda) > 0$, $\sigma_1(\widehat \lambda) = 0$, while for any $\lambda\in(\widehat \lambda, k)$, $\sigma_1(\lambda) < 0$.

Then, if either $k\leq1/(1+\gamma)$ or $k > 1/(1+\gamma)$ but $\lambda\in(\widehat \lambda, k)$, we have $\sigma_{1}(\lambda) < 0$, indicating that $T_n(\lambda) < 0$. Since $\sigma_2(\lambda) > 0$ for $\lambda\in(0, k)$, we have $D_n(\lambda) > 0$. Hence, $(\lambda, v_ \lambda)$ is locally asymptotically stable in system (2.1).

Now, we mainly consider the case when $k > 1/(1+\gamma)$ but $\lambda\in(0, \widehat \lambda)$, in which case $\sigma_1(\lambda) > 0$.

Then, zero is an eigenvalue of $L_n$ if and only if

Defining $d_2 = d_2^{(n)}$, we have ${\rm ker} L = {\rm span}\{\Phi\}$, with $\Phi = (1, \kappa_n)^T\phi_n$, where $\phi_n$ is the eigenfunction $-\Delta$ corresponding to $\tau_n$, and $\kappa_n = (\sigma_1(\lambda)-d_1\tau_n)/\theta > 0$.

Consider the adjoint operator of $L$:

Then, ${\ker} L^* = {\rm span}\{\Phi^*\}$, with $\Phi = (1, \kappa_n^*)^T\phi_n$, where $\kappa_n^* = (d_1\tau_n-\sigma_1(\lambda))/\sigma_2(\lambda) < 0$. Since ${\rm range} \; L = ({\ker} L^*)^{\perp}$, the codimension of ${\rm range} L$ is the same as ${\rm dim} \; L^* = 1$. Thus, $\ker\mathcal F_{(u, v)}(d_2^{(n)}, \lambda, v_ \lambda)$ and $Y/{\rm range} \; \mathcal F_{(u, v)}(d_2^{(n)}, \lambda, v_ \lambda)$ are one-dimensional.

Finally, we show that $\mathcal F_{(d_2, (u, v))}(d_2^{(n)}, \lambda, v_ \lambda)\Phi\notin{\rm range} \; \mathcal F_{(u, v)}(d_2^{(n)}, \lambda, v_ \lambda)$. In fact, we have

and $\big(\mathcal F_{(d_2, (u, v))}(d_2^{(n)}, \lambda, v_ \lambda)\Phi, \Phi^*\big)_Y = (-\tau_n\kappa_n\phi_n, k_n^*\phi_n)_{L^2} = -\tau_n\kappa_nk_n^* > 0$. Thus, we can conclude that $\mathcal F_{(d_2, (u, v))}(d_2^{(n)}, \lambda, v_ \lambda)\Phi\notin{\rm range} \; \mathcal F_{(u, v)}(d_2^{(n)}, \lambda, v_ \lambda)$.

So far, by Theorem 1.7 of ^[8], we have the following results on the local existence of the steady state bifurcation branches:

Theorem 2.1. Suppose that $\frac{1}{4} < \gamma\leq\dfrac{1}{2}$, $k > 1/(1+\gamma)$, $\frac{\theta}{m}\in\bigg(0, \frac{k}{1+k}\bigg(1- \frac{\gamma k}{1+k}\bigg)\bigg)$, $\lambda\in(0, \widehat \lambda)$ and $n$ is a positive integer such that $\tau_n < \sigma_1(\lambda)$ and $d_2^{(n)}\neq d_2^{(m)}$ for any integers $m\neq n$. Then, $(d_2^{(n)}, \lambda, v_ \lambda)$ is a steady state bifurcation point. More precisely, there exists a one-parameter family of non-constant positive solutions $(d_2(s), u(s)(x), v(s)(x))$ of system (2.1) for $|s|$ sufficiently small, where $d_2(s)$, $u(s)$, $v(s)$ are continuous functions, $d_2(0) = d_2^{(n)}$, and $u(s) = \lambda+s\phi_n+o(s)$, $v(s) = v_ \lambda+s\kappa_n\phi_n+o(s)$, where $\kappa_n = (\sigma_1(\lambda)-d_1\tau_n)/\theta > 0$. The zero set of $\mathcal F$ consists of two curves $(d_2(s), u(s)(x), v(s)(x))$ and $(d_2, \lambda, v_ \lambda)$ in a neighborhood of the bifurcation point $(d_2^{(n)}, \lambda, v_ \lambda)$.

Definition 2.2. (^[13]) A component $\mathcal A$ of the set of the non-constant solution of system (2.1) is said to exist globally with respect to $d_2$ if and only if $\overline{\text{Proj}(\mathcal A)}$ contains 0, where $\text{Proj}$ stands for the projection operator from $(0, +\infty)\times X$ to $d_2$ space and the upper bar represents the closure operator in ${\mathbf R}$.

Next, we shall consider the existence of global steady state bifurcation branches.

Theorem 2.3. Let $\mathcal C_1$ be the closure of the set of the non-constant solution of system (2.1) in $(0, +\infty)\times X$. Denote by $\mathcal C_{d_2^{(n)}}$ the component in $\mathcal C_1$ to which $(d_2^{(n)}, \lambda, v_ \lambda)$ belongs. Under the same assumptions of Theorem 2.1, the component $\mathcal C_{d_2^{(n)}}$ exists globally with respect to $d_1$.

Proof. 1) We first argue that, if $(u(x), v(x))$ is a non-negative solution of (2.1), then, either $(u, v)$ is one of $(0, 0)$ and $(k, 0)$, or for $x\in \overline\Omega$, $(u(x), v(x))$ is a non-constant positive solution such that $0 < u(x) < k$ and $0 < v(x) < k(d_2 +\theta d_1)/(\theta d_2)$. This can be obtained by slightly modifying the proof of Lemma 3.5 in ^[12].

2) We are going to use Theorems 2.2 and 2.3 of ^[13] to prove our desired results. Before going, we would like to comment that our results (Theorem 2.3) are similar to but different from Theorem 2.2 of ^[13] since Theorems 2.2 and 2.3 of ^[13] chooses $d_1$ as main bifurcation parameter (fixing $d_2$), while ours chooses $d_2$ as main bifurcation parameter (fixing $d_1$).

Following ^[13], we define the hyperbolic curves $\{C_n\}_{n = 1}^\infty$ in $(0, +\infty)\times(0, +\infty)$ by

Let $\mathcal C_1$ be the closure of the set of the non-constant solution of system (2.1) in $(0, +\infty)\times X$. Denote by $\mathcal C_{d_2^{(n)}}$ the component in $\mathcal C_1$ to which $(d_2^{(n)}, \lambda, v_ \lambda)$ belongs.

To study the existence of global steady state bifurcations, it means to consider how the local bifurcating branch behaves when $d_2$ leaves the critical value $d_2^{(n)}$. By the Rabinowitz's standard global bifurcation theorem ^[9], the component in $(0, +\infty)\times X$ which contains the local bifurcating branch exists globally in the sense that either there is no closed bounded set in the interior of $(0, +\infty)\times X$ which contains the component, or else such a set must contain the bifurcation points other than $(d_2^{(n)}, \lambda, v_ \lambda)$. We shall rule out the possibility of the latter case.

For each fixed $d_1$, we denote the countable set $\{d_2^{(n)}: (d_1, d_2^{(n)})\in\cup_{n = 1}^\infty C_n\backslash\{(d_1, d_2)\in C_n\cap C_m: m\neq n\}\}$ by $\mathcal B_{d_1}$.

Since the non-constant positive solution of system (2.1) is bounded as shown in part 1, by the Rabinowitz's standard global bifurcation theorem, to prove the desired results, it sufficient to rule out the possibility that $\mathcal C_{d_2^{(n)}}$ contains a finite subset $P = \{d_2^{(n)}: (d_2^{(n)}, \lambda, v_ \lambda)\in \mathcal C_{d_2^{(n)}}\}$ of $\mathcal B_{d_1}$. Suppose otherwise. Following ^[13], define

and denote the value of $d_2^{(q)}\in P$ which attains $\overline q$ by $d_2^{(\overline q)}$. Under the assumptions of Theorem 2.1, $d_2^{(n)}\neq d_2^{(m)}$ for any integers $m\neq n$. Thus, the bifurcation point $(d_1, d_2^{(\overline q)})\in C_{\overline q}$ is simple. By Remark 1.2 of ^[13], it follows that if $(d_1, d_2^{(\overline q)})\in C_{\overline q}$, then $\bigg(\overline q^2d_1, \overline q^2d_2^{(\overline q)}\bigg)\in C_1$. This implies that $\bigg(\overline q^2d_1, \overline q^2d_2^{(\overline q)}\bigg)\in C_1$ is simple. In particular, by (2.3) of ^[13], we have

where $T^{\overline q}$ (with $\overline q$ any of the positive integer) is defined in the following way (See also (2.2) of ^[13]):

where $U_0(x)$ is any smooth solution of (2.1) with $(d_1, d_2) = (d_1, d_2^{(n)})$.

Suppose that $\mathcal C_{{(\overline q^2d_2^{(\overline q)})}^{(1)}}$ exists globally (does not contain a bifurcation point different from $(d_2^{(n)}, \lambda, v_ \lambda)$), then $T^{\overline q}\bigg(\mathcal C_{(\overline q^2d_2^{(\overline q))^{(1)}}}\bigg)$ also exists globally (does not contain a bifurcation point different from $(d_1, d_2^{(n)})$, which is impossible since by (2.6), $\mathcal C_{d_2^{(n)}}$ contains a finite subset $P = \{d_2^{(n)}: (d_2^{(n)}, \lambda, v_ \lambda)\in \mathcal C_{d_2^{(n)}}\}$ of $\mathcal B_{d_1}$. A contradiction!

Suppose that $\mathcal C_{{(\overline q^2d_2^{(\overline q)})}^{(1)}}$ does not exist globally, then it must contain a bifurcation point $(d_1, \overline q^2d_2^{(s)})$ for some interger $s\geq2$. However, this is also impossible since it contradicts the fact that $\overline q$ is the maximum value of $q$, since by $T^{\overline q}\bigg(\mathcal C_{{(\overline q^2d_2^{(\overline q)})}^{(1)}}\bigg)$, we can obtain that $\mathcal C_{d_2^{(n)}}$ must contain $(d_2^{(\overline q)s}, \lambda, v_ \lambda)$. Again, the contradiction occurs. Thus, the proof is completed.

3.   The spatial patterns of the shadow system: the time-mapping methods

In this section, we shall consider the existence of non-constant positive steady state solution of the shadow system of the original system, say letting $d_2\rightarrow \infty$.

Rewrite the second equation of system (2.1) in the following way:

Since $0 < u(x) < k$ and $0 < v(x) < k(d_2 +\theta d_1)/(\theta d_2)$, then for fixed $d_1 > 0$, as $d_2\rightarrow \infty$, we have $0\leq u(x)\leq k, 0\leq v(x)\leq k/\theta$. Then, as $d_2\rightarrow \infty$, we have, $\Delta v\rightarrow0$. Therefore, $v\rightarrow \xi$, where $\xi$ is a constant satisfying $0\leq\xi\leq k/\theta$.

Thus, as $d_2\rightarrow \infty$, system (2.1) is reduced to the following limit system:

with $0\leq u(x)\leq k, 0\leq \xi\leq k/\theta$, $\gamma\in(1/4, 1)$. System (3.2) is called shadow system of (2.1).

As in Theorem 2.1, we assume that $1/4 < \gamma\leq1/2$ and $\frac{\theta}{m}\in\bigg(0, \frac{k}{1+k}\big(1- \frac{\gamma k}{1+k}\big)\bigg)$ hold so that system (3.2) has a unique positive constant solution, denoted by $(\lambda, v_ \lambda)$.

We claim that the set of increasing solutions will be able to characterizes all the solutions of (3.2). In fact,

1) The set of increasing solutions can describe all the decreasing solutions of (3.2). Indeed, suppose that system (3.2) has an increasing solution $u = u(x)$ in $(0, \ell)$, then we can check that $u_-(x): = u(\ell-x)$ must be a decreasing solution of system (3.2); Thus, by using increasing solutions, we can construct decreasing solutions;

2) The set of increasing solutions can describe all the solutions of (3.2) (not necessarily limited to the decreasing solutions). Indeed, suppose that system (3.2) has an increasing solution $u = u(x)$ in $(0, \ell)$. Following Lemma 2.1 of ^[13], we define

for $i\ell/n\leq x\leq (i+1)\ell/n$ and $0\leq i\leq n-1$, with $n > 2$; In fact, $u_n(x)$ is constructed by means of successive $(n-1)$-times reflections of $u(nx)\; (0\leq x\leq\ell/n)$ at $i\ell/n\; (i = 1, 2, \cdots, n-1)$. Since (3.2) is an autonomous system and $u_n(x)$ is matched at $i\ell/n\; (i = 1, 2, \cdots, n-1)$. Thus, one can check that $u_n(x)$ is the solution of $d_1/n^2\Delta u+u\big(1-\frac{u}{k}\big)-\frac{mu}{1+u}(1-\frac{\gamma u}{1+u})\xi = 0$. Since the purpose of the paper is to find the non-constant solution of $d_1/\Delta u+u\big(1-\frac{u}{k}\big)-\frac{mu}{1+u}(1-\frac{\gamma u}{1+u})\xi = 0$ for small $d_1$ (See the conditions $f'(\lambda) > \frac{d_1\pi^2}{\ell^2}$ in Theorem 3.3). Thus, if we are able to prove the existence of the increasing solution of (3.2) in $(0, \ell)$, we can construct other non-monotone solutions $u_n(x)$ by using $d_1/n^2$ to replace $d_1$;

In the following, we shall only concentrate on the increasing solutions of system (3.2).

We define

and introduce the energy functional $I(x): = \frac{d_1}{2}(u_x(x))^2+F(u(x))$. Then, $I'(x) = d_1u''(x)+f(u)$ in $(0, \ell)$. If $u = u(x)$ is a non-constant increasing solution of (3.2), then $I'(x)\equiv0$ in $(0, \ell)$ and $F(u(x)) < F(u(0)) = F(u(\ell))$, since $u'(0) = u'(\ell) = 0$. In particular, $F(u)$ must attain its local minimal value at a point in $(u(0), u(\ell))$.

We rewrite $f(u)$ as

Clearly, $p(0) = 1/m$, $p(+\infty) = -\infty$. We have the following results on the function $p(u)$:

Lemma 3.1. 1) Suppose that $k\leq1/(1+\gamma)$ holds. Then, $p(u)$ is strictly decreasing in $(0, \infty)$ and there exists a unique $u^*\in(0, \infty)$, such that for any $u\in(0, u^*)$, $p(u) > 0$ for any $u\in(0, u_*)$, while $p(u) < 0$ for any $u\in(u^*, \infty)$. In particular, $p(u)$ attains its local maximal value at $u = 0$.

2) Suppose that $k > 1/(1+\gamma)$ holds. Then, there exists $u_*, u^*\in(0, \infty)$, with $u_* < u^*$, such that for any $u\in(0, u_*)$, $p'(u) > 0$ for any $u\in(0, u_*)$, while $p'(u) < 0$ for any $u\in(u_*, \infty)$. In particular, $p(u)$ attains its maximal value at $u = u_*$. Moreover, for any $u\in(0, u^*)$, $p(u) > 0$ for any $u\in(0, u_*)$, while $p(u) < 0$ for any $u\in(u^*, \infty)$.

Proof. Note that

Analyzing the properties of $p'(u)$, we can obtain the desired results.

From (3.4) and Lemma 3.1, it follows that: if either $k\leq1/(1+\gamma)$ or $k > 1/(1+\gamma)$ but $\xi\geq p(u_*)$ holds, $F(u)$ will never its minimal value point in $(0, \infty)$. Therefore, system (3.2) will never have positive non-constant solutions. Thus, to expect non-constant positive solutions, we need to assume that $k > 1/(1+\gamma)$ and $\xi\in\big(1/m, p(u_*)\big)$. In this case, by (3.4), $f(u) = 0$ have two roots in $(0, \infty)$, denoted by $u_{1, 2}(\xi)$, satisfying $0 < u_1(\xi) < u_* < u_2(\xi) < u^*$; And $f(u) < 0$ for $0 < u < u_1(\xi)$ and $u_2(\xi) < u < +\infty$, while $f(u) > 0$ for $u_1(\xi) < u < u_2(\xi)$. Therefore, $F(u)$ is decreasing in $(0, u_1(\xi))\cup(u_2(\xi), \infty)$, and increasing in $(u_1(\xi), u_2(\xi))$. In particular, $F(u)$ takes its local minimum value at $u = u_1(\xi)$. Define

where $\kappa(\xi)\in(0, u_1(\xi))$, such that $F(\kappa(\xi)) = F(u_2(\xi))$. Then, for any $u(0)\in(\kappa_0(\xi), u_1(\xi))$, there exists a unique $u^0(\xi)\in(u_1(\xi), u_2(\xi))$, such that $F(u(0)) = F(u^0(\xi))$.

Since $I'(x)\equiv0$, we have $I(x)\equiv F(\alpha)$, where we write $u(0) = \alpha$ and $u(\ell) = \beta$. Since $u'(x)$ is assumed to be positive, we have

Integrating (3.5) from $0$ to $\ell$, by $u(0) = \alpha$ and $u(\ell) = \beta$, we have,

Regarding $\ell$ in (3) as a function of $\alpha$, we have the following results:

Lemma 3.2. Suppose that $\dfrac{1}{4} < \gamma\leq\dfrac{1}{2}$, $\xi\in(1/m, p(u_*))$ and $k > \dfrac{1}{1+\gamma}$ hold. Then, for any $\alpha\in(\kappa_0(\xi), u_1(\xi))$, $\dfrac{d\ell}{d \alpha} < 0$.

Proof. Step 1: Following ^[10], for a given number $u\in(\alpha, \beta)$, define $u = \rho(s)$ by the relation

Differentiating the first equation of (3.6) with respect to $s$, we have $f(\rho(s))\rho'(s) = s$. By the second equation of (3.6), we have $\rho(s) > u_1(\xi)$ for $s > 0$, $\rho(0) = 0$ and $\rho(s) < u_1(\xi)$ for $s < 0$. This together with $f(\rho(s))\rho'(s) = s$, implies that $\rho'(s) > 0$ for all $s\neq0$, where we use the facts that $f(u) > 0$ for $u\in(u_1(\xi), u_2(\xi))$ and $f(u) < 0$ for $u\in(\kappa_0(\xi), u_1(\xi))$. Then, $s = \rho^{-1}(u)$ is well defined and is strictly increasing in $(\alpha, \beta)$. By ^[10] (pages of 314–315), we have

and $\rho'''(s) = - \frac{\rho'(s)}{f^4(u)}z(u)$, where

In particular, $\rho'(0) = \frac{1}{\sqrt{f'(u_1(\xi))}}$ and

Step 2: Let $\omega > 0$ be given by

which implies from (3.6) that $2(F(\alpha)-F(u)) = \omega^2-s^2$ and

where the last equality was obtained by making the change of variable $s = -\omega\cos t$, $0\leq t\leq\pi$. Regarding $\omega$ in (3.10) as a function of $\alpha$ in $(\kappa_0(\xi), u_1(\xi))$, by $f(\alpha) < 0$ and $\omega(\alpha) > 0$, we have $\omega'(\alpha) = f(\alpha)/\omega(\alpha) < 0$. Thus, the sign of $\frac{d\ell}{d \alpha}$ is opposite to the sign of $\frac{d\ell}{d\omega}$. By (3.11), we have

After elementary calculations, we have

Since $m\xi > 1$ and $k > 1/(1+\gamma)$, we have $f''(0) = 2m\xi(\gamma+1)- \frac{2}{k} > 0$. Note that $\gamma\leq1/2$ holds. Then, $f'''(u) < 0$ for all $u > 0$. Thus, $f''(u)$ is decreasing in $(0, \infty)$. Since $f''(0) > 0$, by $f''(\infty) = -2/k$ and the properties of $f'''(u)$, $f''(u) = 0$ has a unique root $\mu_*$ in $(0, \infty)$, such that $f''(u) > 0$ in $(0, \mu_*)$ and $f''(u) < 0$ in $(\mu_*, \infty)$. We argue that $f'(\mu_*) > 0$. Otherwise, since $f'(u)$ attains its maximum value at $u = \mu_*$, we have $f'(u)\leq0$ for all $u > 0$. Since $f'(0) = 1-m\xi < 0$, we have $f'(u) < 0$. Since $f(0) = 0$, thus $f(u) < 0$, which is again impossible. Thus, $f'(\mu_*) > 0$. By $f'(0) < 0$ and $f'(\infty) = -\infty$, $f'(u) = 0$ has two roots $\mu_{1}$ and $\mu_2$ in $(0, \infty)$ with $\mu_1 < \mu_* < \mu_2$. In particular, $0 < \mu_1 < u_1(\xi) < \mu_2 < u_2(\xi)$.

Step 3: Since $f(u_1(\xi)) = 0$, we have $z(u_1(\xi)) = 0$. For $u\in(0, \mu_1]$, we have $f(u) < 0, \; f'(u)\leq0, f''(u) > 0, F(u)-F(u_1(\xi)) > 0$. Then, for $u\in(0, \mu_1]$, we have $z(u) < 0$. For $u\in[\mu_2, u_2(\xi))$, we have $f(u) > 0, \; f'(u)\leq0, f''(u) < 0, F(u)-F(u_1(\xi)) > 0$. Then, for $u\in[\mu_1, u_2(\xi))$, we have $z(u) < 0$. Next, we show that for $u\in(\mu_1, \mu_2)$, $z(u) < 0$. A straightforward calculation shows that,

where $z_1(u): = 5(f''(u))^2-3f'(u)f'''(u)$. For $u\in(\mu_1, \mu_2)$, we have $f'(u) > 0$. By $f'''(u) < 0$, we have $z_1(u) > 0$. On the other hand, by (3.9) and (3.13), we have $\rho'''(0) = z_1(u_1(\xi))/(12\sqrt{f'(u_1(\xi))^7}) > 0$. This implies that for $s$ close to zero, $\rho'''(s) > 0$. Hence, by $\rho'''(s) = -\rho'(s)z(u)/f^4(u)$, for $u$ close to $u_1(\xi)$, but not equal to $u_1(\xi)$, we have $z(u) < 0$. As in ^[10], suppose for contradiction that there exists a $\widehat{u}\in(u_1(\xi), \mu_2)$ such that $z(\widehat{u}) = 0$ and $z(u) < 0$ for $u\in(u_1(\xi), \widehat{u})$, then $z'(\widehat{u})\geq0$ and $f'(\widehat{u}) > 0$ and $z_1(\widehat{u}) > 0$, which contradicts (3.13). Thus, $z(u) < 0$ for all $u\in(u_1(\xi), u_2(\xi))$. Similarly, we can show that $z(u) < 0$ for all $u\in(0, u_1(\xi))$. Finally, by $\rho'''(s) = -\rho'(s)z(u)/f^4(u)$ and (3.12), for $\alpha\in(\kappa_0(\xi), u_1(\xi))$, $d\ell/d \alpha < 0$.

So far, we are in the position to state our results on the existence of positive non-constant solutions of the shadow system:

Theorem 3.3. Suppose that $\dfrac{1}{4} < \gamma\leq\dfrac{1}{2}$, $f'(\lambda) > \frac{d_1\pi^2}{\ell^2}$ and $k > \dfrac{1}{1+\gamma}$ hold, where $\lambda$ is the first component of $(\lambda, v_ \lambda)$. Then, the shadow system (3.2) admits at least one strictly increasing solution, and the same number of strictly decreasing solutions. The corresponding $\xi$ satisfies $\xi\in(\xi_-, \xi_+)$, for some $\xi_-, \xi_+\in(0, p(u_*))$.

Proof. For $\xi\in\big(1/m, p(u_*)\big)$, since $\rho'(0) = \sqrt{d_1}\pi/f'(u_1(\xi))$, we have $\lim_{ \alpha\rightarrow u_1(\xi)}\ell(\alpha) = \rho'(0) = \sqrt{d_1}\pi/\sqrt{f'(u_1(\xi))} = :\ell_0$. Then, for $\ell < \ell_0$, $d_1\Delta u+f(u) = 0$ has no non-constant solutions, while for $\ell > \ell_0$ it does have. On the other hand, from the argument in the proof of Theorem 5 in ^[10], it follows that $\lim_{ \alpha\rightarrow \kappa_0(\xi)}\ell(\alpha) = +\infty$ and $u$ is non-degenerate.

As in the Step 3 of the proof of Lemma 3.2, we have shown that $f'(u_1(\xi))$ attains its unique maximal value (positive) at $u_1(\xi) = \mu_*$ for some $\mu_*\in(\mu_1, \mu_2)$. Define $\ell_* = \sqrt{d_1}\pi/\sqrt{f'(\mu_*)}$. Then, $f'(\mu_*) = d_1\pi^2/\ell_*^2$ and for any $\ell > \ell_*$, there have two numbers $u_{\ell_-} < u_{\ell_+}$, such that $\delta_*\in(u_{\ell_-}, u_{\ell_+})\subset\big(0, u_*\big)$ and $f'(u_{\ell_-}) = f'(u_{\ell_+}) = d_1\pi^2/\ell^2$. Let $\xi_-: = p(u_{\ell_-})$ and $\xi_+: = p(u_{\ell_+})$. Then, for $\xi\in(\xi_-, \xi_+)$, we have $f'(u_1(\xi)) > d_1\pi^2/\ell^2$; Thus, there must be a unique, increasing and non-degenerate solution $u = u_\ell(x)$ to $d_1\Delta u+f(u) = 0$ subject to Neumann boundary conditions. Recall that under our assumptions, $u_1(\xi) = \lambda$ if and only if $\xi = v_ \lambda$. Then, if $f'(\lambda) > d_1\pi^2/\ell^2$, we have $f'(\mu_*)\geq f'(\lambda) > d_1\pi^2/\ell^2$, which implies that $\ell > \ell_*$ and $\xi\in(\xi_-, \xi_+)$. In particular, $\lambda\in(u_{\ell_-}, u_{\ell_+})$. As $\ell\rightarrow \ell_*$, the set $(\xi_-, \xi_+)$ shrinks to an empty set, and as $\ell\rightarrow \infty$, $(\xi_-, \xi_+)$ expands to $(1/m, p(u_*))$. For fixed $\ell > \ell_*$ and a given $\xi\in(\xi_-, \xi_+)$, we write

It follows from the non-degeneracy of $u_\ell$ that both $\alpha_\ell(\xi)$ and $\mathcal M(\xi)$ are continuous functions of $\xi$. If $\xi\in(\xi_-, \xi_+)$ is close to $\xi_-$, then $u_1(\xi)$ is close to $u_{\ell_-}$. Since $\ell = \sqrt{d_1}\pi/\sqrt{f'(u_{\ell_-})}$, we know that $\alpha_\ell$ is close to $u_1(\xi)$ and then close to $u_{\ell_-}$ too. Furthermore, the solution $u_\ell(x)$ is nearly a constant. Then, by $u_{\ell_-} < \lambda$, $\mathcal M(\xi)$ is close to $u_{\ell_-}$ and then we have $\mathcal M(\xi) < \mathcal M(\lambda) = \theta/m$, where the integrant function in $\mathcal M(\xi)$ is increasing with respect to $u$ when $\gamma\leq1/2$. Similarly, for those $\xi$ close to $u_{\ell_+}$, by $u_{\ell_+} > \lambda$, we have $\mathcal M(\xi) > \mathcal M(\lambda) = \theta/m$. Then, by the continuity of $\mathcal M(\xi)$, there must have some $\xi\in(\xi_-, \xi_+)$ such that $\mathcal M(\xi) = \theta/m$.

4.   Conclusions

In this paper, we are mainly concerned with the existence of the spatial patterns of the diffusive system with plant-herbivore interactions with toxin-determined functional responses. By using the standard steady state bifurcation theorem, we are able to show that, for any fixed $d_1 > 0$, if the diffusion rate $d_2$ is larger than certain critical value, then the original non-degenerate system (both $d_1$ and $d_2$ are positive and finite) will have non-constant positive steady state solutions. Then, we continue to consider the limiting case when $d_2$ tends to infinity. The corresponding limiting system is called the shadow system of the original system. By using the time-mapping methods, we have showed that the shadow system has non-constant positive steady state solutions under certain conditions. In the future, we shall concentrate ourselves on the dynamics of the predator-prey system with discontinuous reaction terms, which tends to be much more interesting than the current counterpart (see for example, ^[14,15]).

Acknowledgments

The author would like to thank the reviewers for their kind suggestions and comments which significantly help improve the quality of the paper.

Conflict of interest

The authors declare there is no conflicts of interest.