Global Dynamics and Existence of Traveling Wave Solutions for A Three-Species Models

Fanfan Li School of Mathematical Sciences, University of Jinan, Jinan, Shandong 250022, P R China. Zhenlai Han School of Mathematical Sciences, University of Jinan, Jinan, Shandong 250022, P R China. Ting-Hui Yang¹¹1Research was partially supported by Ministry of Science and Technology, Taiwan (R O C). Department of Mathematics, Tamkang University, Tamsui Dist., New Taipei City 25137, Taiwan.

Abstract

In this work, we investigate the system of three species ecological model involving one predator-prey subsystem coupling with a generalist predator with negative effect on the prey. Without diffusive terms, all global dynamics of its corresponding reaction equations are proved analytically for all classified parameters. With diffusive terms, the transitions of different spatial homogeneous solutions, the traveling wave solutions, are showed by higher dimensional shooting method, the Wazewski method. Some interesting numerical simulations are performed, and biological implications are given.

2010 Mathematics Subject Classification. Primary: 37N25, 35Q92, 92D25, 92D40.
Keywords : Two predators-one prey system, extinction, coexistence, global asymptotically stability, traveling wave solutions, Wazewski principle.

1 Introduction

In this work, we consider an ecological system of three species with diffusion as follows,

\left\{\begin{aligned} \partial_{t}u&=r_{1}u(1-u)-a_{12}uv-a_{13}uw,\\ \partial_{t}v&=r_{2}v(1-v)+a_{21}uv,\\ \partial_{t}w&=d\Delta w-\mu w+a_{31}uw,\end{aligned}\right.

(1.1)

where parameters $d$ is the diffusive coefficient for species $w$ , $r_{1}$ and $r_{2}$ are the intrinsic growth rates of species $u$ and $v$ respectively, and $\mu$ is the death rate of the predator $w$ . The nonlinear interactions between species is the Lotka-Volterra type interactions between species where $a_{ij}(i<j)$ is the rate of consumption and $a_{ij}(i>j)$ measures the contribution of the victim (resource or prey) to the growth of the consumer [19]. Here, the species $u$ is the nutrient resource of the predator-prey system, the species $v$ is called the generalist predator which can take advantage of various resources from two trophic levels, and the species $w$ is called the specialist predator which has a limited diet from $u$ . To simplify the analysis, we only consider the species $w$ has diffusion effect.

System (1.1) without diffusive terms is given by the following system of three ODEs:

\left\{\begin{aligned} \dot{u}&=r_{1}u(1-u)-a_{12}uv-a_{13}uw,\\ \dot{v}&=r_{2}v(1-v)+a_{21}uv,\\ \dot{w}&=-\mu w+a_{31}uw.\end{aligned}\right.

(1.2)

The whole system can be seen as a predator-prey subsystem, $u$ - $w$ subsystem, coupled with an extra species $v$ with negative effect on species $u$ . The model can be used to describe parasitism, consumption or predation in the community of plants species [7]. It is well-known that the ecological principle of competitive exclusion holds for the following classical two predators-one prey model [17],

\left\{\begin{aligned} \dot{u}&=r_{1}u(1-u)-a_{12}uv-a_{13}uw,\\ \dot{v}&=-\mu_{1}v+a_{21}uv,\\ \dot{w}&=-\mu_{2}w+a_{31}uw.\end{aligned}\right.

However, by comparing these two models, system (1.2) is a modified two-predators one-prey model where one, $v$ , is a generalist predator and another one, $w$ , is a specialist predator. There are a fundamental difference, in Section 2, where we show the positive equilibrium of system (1.2) can not only exist but it is also globally asymptotically stable, that is, two predators can co-exist.

Reaction-diffusion systems are often characterized by the existence of spatial homogeneous equilibria when the diffusion terms vanish. If there are more than one equilibrium, then we can expect a possible transition between them. These transitions are described by reaction-diffusion waves. Propagations of flames, migration of biological species, or tumor growth are among many examples of such phenomena [18, 21]. In the PDE perspective, the existence of traveling wave solutions for reaction-diffusion systems in an important and interesting subject which has attracted considerable attentions [2, 3, 4, 8, 12, 9, 10, 11, 22]. The phenomena of traveling wave solutions of reaction-diffusion systems have been widely studied [21] from the single equation with nonlinearity in monostable type [6] or bistable type [13] to monotone systems [14]. There have been great successes in the existence, uniqueness, stability and spreading speed of traveling wave solutions of monotone system [1, 14, 20].

Unfortunately, our system which has an important nonlinear interaction, predator-prey type, between different species is non-monotone. In the past three decades, by using different methods including the shooting method, Conley index and upper-lower solutions method, the existence of traveling wave solutions has been established for various predator-prey systems. See Dunbar [2, 3, 4], Gardner and Smoller [8], Jones et al. [12], Huang [10, 11], Hsu et al. [9], Lin et al. [15], and references cited therein. In this work, we will use the so-called higher dimensional shooting method, Wazewski method to show the existence of positive traveling wave solutions from one unstable equilibrium to a stable one. Here we briefly describe this framework of the shooting method.

To show the existence of traveling wave solution, by using the moving coordinates, the reaction-diffusion system is transformed into a ODE system, and the existence of traveling wave solutions connecting two different equilibria is equivalent to a heteroclinic orbit of the corresponding ODE system. we analyze the structure of unstable manifold of the unstable equilibrium first. Then we construct a variant of Wazewski set $\Sigma$ with the unstable equilibrium as its boundary point and also containing the stable equilibrium. Third, dynamics of the system on all boundary of $\Sigma$ should be clarified. Next, pick up a curve contained in the unstable manifold with two end points on the “exit set” of boundary of $\Sigma$ . It is clear that all solutions with initial conditions on this curve will attend to the unstable equilibrium as $t\to-\infty$ . Then show that there exists a particular point on this curve, and the solution starting from this point will stay in the interior of $\Sigma$ for all $t\geq 0$ . Finally, we define a nonempty subset of $\Sigma$ which contains the point stayed in $\Sigma$ for all positive time under the action of the ODE system. Then we can get our main result by constructing a Lyapunov function and use the LaSalle’s invariance principle on this non-empty set.

In this work, our main contributions are as follows. First, for the corresponding reaction system (1.2) we clarify completely the existence, non-existence, and all asymptotically states and their global stabilities are investigated theoretically. Secondly, we show the existence of traveling wave solutions are obtained for a particular three species ecosystem with predator-prey interaction. We use the Wazewski principle to show the existence of traveling wave solutions. Thought, the method is similar to [10, 11], our system is three dimensional. Third, numerical simulations are performed for some interesting initial functions. Finally, some biological interpretations are given.

The rest of this article is organized in the following manner. In Section 2, we first consider the corresponding reaction equations of (1.1) which is a system of three ODEs. The existence of boundary equilibria and coexistence equilibrium are obtained with some conditions. Moreover, we find the necessary and sufficient condition of global asymptotic stability of the positive equilibrium. In Section 3, by using shooting method, the existence and nonexistence of traveling wave solutions of (1.1) are obtained. In Section 4, the numerical simulations are performed and presented. Finally, some remarks and discusses in biological meanings are also given in the last section.

2 Global Dynamics of the Corresponding ODE System

In this section, we investigate dynamics of the ODE system (1.2) and the essential assumptions to guarantee the existence and local stabilities of all equilibria. Moreover, two extinction results and the global stability of positive equilibrium are showed.

2.1 Preliminaries

It is easy to see that $uv$ -, $uw$ -, and $vw$ -planes are invariant subspaces of (1.2). Hence solutions of (1.2) will be positive/non-negative if they start from a positive/non-negative point. Moreover, we can show that solutions are bounded.

Lemma 2.1.

The solutions of system (1.2) are bounded.

Proof.

From the first equation in system (1.2) we have

\frac{du}{dt}=r_{1}u(1-u)-a_{12}uv-a_{13}uw\leqslant r_{1}u(1-u),

so that the comparison principle implies that

\limsup_{t\rightarrow\infty}u\leqslant 1,

i.e. $0<u\leqslant 1$ . Similarly we can get $0<v\leqslant 1+\frac{a_{21}}{r_{2}}$ .
Let $M=r_{1}+\mu$ and $D=\mu$ . From the first and third equation in (1.2) we have

	$\displaystyle\frac{d}{dt}\left(u+\frac{a_{13}}{a_{31}}w\right)$	$\displaystyle=r_{1}u(1-u)-a_{12}uv-a_{13}uw-\frac{a_{13}}{a_{31}}w+a_{13}uw$
		$\displaystyle\leq r_{1}u-\frac{a_{13}}{a_{31}}w$
		$\displaystyle=Mu-D(u+\frac{a_{13}}{a_{31}}w)$
		$\displaystyle\leq M-D(u+\frac{a_{13}}{a_{31}}w).$

Using the comparison principle again, we have

\limsup_{t\rightarrow\infty}\left(u+\frac{a_{13}}{a_{31}}w\right)\leq\frac{M}{D}.

∎

2.2 Assumptions and Two Extinction Results

From now on, we always make the assumptions,

(H1)

$r_{1}>a_{12}$ ,
(H2)

$a_{31}>\mu$ ,

which will be used in the rest of the article, because of the following two extinction results.

Lemma 2.2.

If $r_{1}\leq a_{12}$ , then $\lim_{t\rightarrow\infty}u(t)=0$ and $\lim_{t\rightarrow\infty}w(t)=0$ .

Proof.

It is easy to see that if $\lim_{t\rightarrow\infty}u(t)=0$ then $\lim_{t\rightarrow\infty}w(t)=0$ , sequently. Hence we only show that the first limit holds. Two cases, $r_{1}<a_{12}$ and $r_{1}=a_{12}$ , are considered. It is easy to see that $v(t)\geq 1$ eventually by comparison principle, since

\dot{v}=r_{2}v(1-v)+a_{21}uv\geq r_{2}v(1-v).

For $r_{1}<a_{12}$ and $t$ large enough, we have

\frac{\dot{u}}{u}=r_{1}(1-u)-a_{12}v-a_{13}w<r_{1}-a_{12}<0.

Hence we have $u(t)<ce^{(r_{1}-a_{12})t}\to 0$ as $t\rightarrow\infty$ .

For $r_{1}=a_{12}$ , we obtain that $u(t)$ is decreasing with respect to $t$ , and claim that $\lim_{t\to\infty}u(t)=0.$ Suppose to the contrary that $\lim_{t\to\infty}u(t)=\xi>0$ . By Markus limiting theorem [16], we have

	$\displaystyle 0=\liminf_{t\to\infty}\dot{u}(t)$	$\displaystyle=\liminf_{t\to\infty}\big{(}r_{1}u(t)(1-u(t))-r_{1}u(t)v(t)-a_{13}u(t)w(t)\big{)}$
		$\displaystyle=r_{1}\xi(1-\xi)-r_{1}\xi\limsup_{t\to\infty}v(t)-a_{13}\xi\limsup_{t\to\infty}w(t)$
		$\displaystyle=\xi\big{(}r_{1}(1-\xi)-r_{1}\limsup_{t\to\infty}v(t)-a_{13}\limsup_{t\to\infty}w(t)\big{)}=0,$

which implies

\displaystyle r_{1}(1-\xi)=r_{1}\limsup_{t\to\infty}v(t)+a_{13}\limsup_{t\to\infty}w(t).

(2.1)

However, by the second equation of (1.2) and the comparison principle, it is clear that $\liminf_{t\to\infty}v(t)\geq 1$ . Hence, by (2.1), we obtain that

a_{13}\limsup_{t\to\infty}w(t)=r_{1}(1-\xi)-r_{1}\limsup_{t\to\infty}v(t)\leq r_{1}(1-\xi)-r_{1}<0

which contradicts to the positivity of $w(t)$ . Hence we complete the proof. ∎

Lemma 2.3.

If $a_{31}\leq\mu$ , then $\lim_{t\rightarrow\infty}w(t)=0$ .

Proof.

Similarly, two cases, $a_{31}<\mu$ and $a_{31}=\mu$ , are considered. For case $a_{31}<\mu$ , we have

\frac{\dot{w}}{w}=a_{31}u-\mu\leq a_{31}-\mu<0

holds. Hence $w(t)\to 0$ as $t\rightarrow\infty$ .

For case $a_{31}=\mu$ ,

\dot{w}=w(a_{31}u-\mu)\leq w(a_{31}-\mu)=0.

Hence $w$ is decreasing, and we claim that $\lim_{t\to\infty}w(t)=0$ . Suppose to the contrary that $\lim_{t\rightarrow\infty}w(t)=\xi>0$ . By comparison principle again, we have

	$\displaystyle 0=\liminf_{t\rightarrow\infty}\dot{w}(t)$	$\displaystyle=\liminf_{t\rightarrow\infty}\big{(}-\mu w(t)+a_{31}u(t)w(t)\big{)}$
		$\displaystyle=a_{31}\xi\liminf_{t\rightarrow\infty}u(t)-\mu\xi$
		$\displaystyle=a_{31}\xi\big{(}\liminf_{t\rightarrow\infty}u(t)-1\big{)}=0,$

which implies $\lim_{t\to\infty}u(t)=1$ and $\lim_{t\to\infty}\dot{u}(t)=0$ . However, by considering the first equation of (1.2), we obatin

0=\lim_{t\rightarrow\infty}\dot{u}(t)=\lim_{t\to\infty}\big{(}r_{1}u(1-u)-a_{12}uv-a_{13}uw\big{)}\leq-a_{13}\xi<0,

where this is a contradiction. Hence we have $\lim_{t\rightarrow\infty}w(t)=0$ , and the proof is complete. ∎

Biologically, these two results can be easily interpreted in the biological point of view. From the first equation of (1.2), species $u$ have two negative effects from $v$ and $w$ , respectively. To sustain the negative effect of species $v$ , $r_{1}>a_{12}$ , is necessary for survival of $u$ and supporting for $w$ in Lemma 2.2. Alternatively, for Lemma 2.3, if the mortality rate $\mu$ of species $w$ is greater than the benefit getting from species $u$ , the conversion rate $a_{31}$ , then $w$ will die out eventually. Whenever $r_{1}\leq a_{12}$ or $\mu\leq a_{31}$ hold, then system (1.2) is reduced to a one- or two-dimensional subsystem of (1.2) which is well studied by classical results like Poincare-Bendixson Theorem. Hence we make assumptions (H1) and (H2).

2.3 Equilibria and Stability in $\mathbb{R}^{3}$

By straightforward calculation, we obtain that there are one trivial equilibrium $E_{0}=(0,0,0)$ , and four semi-trivial equilibria, $E_{1}=(1,0,0)$ , $E_{2}=(0,1,0)$ ,

	$\displaystyle E_{12}$	$\displaystyle=(u^{}_{12},v^{}_{12},0)=\left(\frac{r_{2}(r_{1}-a_{12})}{r_{1}r_{2}+a_{12}a_{21}},\frac{r_{1}r_{2}+r_{1}a_{21}}{r_{1}r_{2}+a_{12}a_{21}},0\right)\quad\text{ and }$
	$\displaystyle E_{13}$	$\displaystyle=(u^{}_{13},0,w^{}_{13})=\left(\frac{\mu}{a_{31}},0,\frac{r_{1}(a_{31}-\mu)}{a_{13}a_{31}}\right),$

of system (1.2). Here $u^{*}_{12}$ and $v^{*}_{12}$ satisfy the equations,

\left\{\begin{aligned} r_{1}(1-u^{*}_{12})-a_{12}v^{*}_{12}&=0,\\ r_{2}(1-v^{*}_{12})+a_{21}u^{*}_{12}&=0,\end{aligned}\right.

(2.2)

and $u^{*}_{13}$ and $w^{*}_{13}$ have the forms,

\displaystyle u^{*}_{13}=\mu/a_{13}\quad\text{ and}\quad r_{1}(1-u^{*}_{13})=a_{13}w^{*}_{13}.

(2.3)

It is obvious that the equilibria, $E_{0}$ , $E_{1}$ and $E_{2}$ , always exist without any restriction. By contrast, the equilibria $E_{12}$ and $E_{13}$ exist if assumptions (H1) and (H2) hold, respectively. The positive equilibrium

	$\displaystyle E_{*}$	$\displaystyle=(u_{},v_{},w_{*})$
		$\displaystyle=\left(\frac{\mu}{a_{31}},1+\frac{a_{21}\mu}{a_{31}r_{2}},\frac{r_{1}r_{2}a_{31}-r_{1}r_{2}\mu-a_{12}a_{31}r_{2}-a_{12}a_{21}\mu}{a_{13}a_{31}r_{2}}\right)$		(2.4)

exists if

(H3)

$r_{1}r_{2}a_{31}-r_{1}r_{2}\mu-a_{12}a_{31}r_{2}-a_{12}a_{21}\mu>0.$

holds. It is hard to see clearly the biological meanings of assumption (H3) because of the complicated form. However, it is easy to see that the inequality $r_{1}\leq a_{12}$ implies that (H3) does not hold, that is, (H3) is a sufficient conditions of (H1). Similarly, assumption (H3) is also a sufficient conditions of (H2). Biologically, there are two key points for the existence of positive equilibrium $E_{*}$ . One is the survival of species $u$ by its $r$ -strategy to overcome the negative effect from species $v$ , and another one is species $w$ should overcome its mortality by getting benefit from species $u$ .

By direct computations, we have the Jacobian matrix of system (1.2) given by

J=\begin{bmatrix}r_{1}-2r_{1}u-a_{12}v-a_{13}w&-a_{12}u&-a_{13}u\\ a_{21}v&r_{2}-2r_{2}v+a_{21}u&0\\ a_{31}w&0&-\mu+a_{31}u\end{bmatrix}.

(2.5)

(i)

It is clear that

$J(E_{0})=\begin{bmatrix}r_{1}&0&0\\ 0&r_{2}&0\\ 0&0&-\mu\end{bmatrix}$

has two positive eigenvalues and one negative eigenvalue, and $E_{0}$ is saddle.
(ii)

Evaluating (2.5) at $E_{1}$ implies that

$J(E_{1})=\begin{bmatrix}-r_{1}&-a_{12}&-a_{13}\\ 0&r_{2}+a_{21}&0\\ 0&0&-\mu+a_{31}\end{bmatrix}$

has two positive eigenvalues and one negative eigenvalue. Similarly, we can obtain the matrix

$J(E_{2})=\begin{bmatrix}r_{1}-a_{12}&0&0\\ a_{21}&-r_{2}&0\\ 0&0&-\mu\end{bmatrix}$

which is stable if assumption (H1) does not hold. Actually, we can show that, by Markus limiting theorem [16], equilibrium $E_{2}$ is globally asymptotically stable if (H1) does not hold.

(iii)

The Jacobian evaluated at $E_{12}$ is given by

J(E_{12})=\begin{bmatrix}-r_{1}u^{*}_{12}&-a_{12}u^{*}_{12}&-a_{13}u^{*}_{12}\\ a_{21}v^{*}_{12}&-r_{2}v^{*}_{12}&0\\ 0&0&-\mu+a_{31}u^{*}_{12}\end{bmatrix}.

(2.6)

It is easy to see that there are two eigenvalues, $\lambda_{1}$ and $\lambda_{2}$ , corresponding to the upper-left $2\times 2$ submatrix of (2.6) and one eigenvalue, $\lambda_{3}=-\mu+a_{31}u^{*}_{12}$ , with

	$\displaystyle\lambda_{1}+\lambda_{2}$	$\displaystyle=-(r_{1}u^{}_{12}+r_{2}v^{}_{12})<0,$
	$\displaystyle\lambda_{1}\lambda_{2}$	$\displaystyle=(r_{1}r_{2}+a_{12}a_{21})u^{}_{12}v^{}_{12}>0.$

Thus the matrix has two negative eigenvalues and one positive eigenvalue if $\mu<a_{31}u^{*}_{12}$ , or three negative eigenvalues if $\mu>a_{31}u^{*}_{12}$ .

(iv)

The Jacobian evaluated at $E^{*}_{13}$ is given by

J(E^{*}_{13})=\begin{bmatrix}-r_{1}u^{*}_{13}&-a_{12}u^{*}_{13}&-a_{13}u^{*}_{13}\\ 0&r_{2}+a_{21}u^{*}_{13}&0\\ a_{31}w^{*}_{13}&0&0\end{bmatrix}.

(2.7)

Observing the form of the Jacobian matrix (2.7), there are one positive eigenvalue, $\lambda_{3}=r_{2}+a_{21}u^{*}_{13}$ , and two eigenvalues, $\lambda_{1}$ and $\lambda_{2}$ , which are obtained by removing the second column and the second row. Although it is obvious that

	$\displaystyle\lambda_{1}+\lambda_{2}$	$\displaystyle=-r_{1}u^{*}_{1}3<0\quad\text{ and}$
	$\displaystyle\lambda_{1}\lambda_{2}$	$\displaystyle=a_{1}3a_{3}1u^{}_{1}3w^{}_{1}3>0,$

equilibrium $E_{13}$ is saddle.

(v)

The Jacobian evaluated at $E_{*}$ is given by

J(E_{*})=\begin{bmatrix}-r_{1}u_{*}&-a_{12}u_{*}&-a_{13}u_{*}\\ a_{21}v_{*}&-r_{2}v_{*}&0\\ a_{31}w_{*}&0&0\end{bmatrix}.

Then the characteristic equation is

	$\displaystyle\lambda^{3}+(r_{1}u_{}+r_{2}v_{})\lambda^{2}+$	$\displaystyle(r_{1}r_{2}u_{}v_{}+a_{12}a_{21}u_{}v_{}+a_{13}a_{31}u_{}w_{})\lambda$
		$\displaystyle+r_{2}a_{13}a_{31}u_{}v_{}w_{*}=0.$

It is obviously that, by Routh-Hurwitz criterion, the real parts of three roots of the characteristic equation are all negative if and only if

(r_{1}u_{*}+r_{2}v_{*})(r_{1}r_{2}u_{*}v_{*}+a_{12}a_{21}u_{*}v_{*}+a_{13}a_{31}u_{*}w_{*})>r_{2}a_{13}a_{31}u_{*}v_{*}w_{*}

which is clearly true. Hence the positive equilibrium $E_{*}$ is stable whenever it exists.

Let us summarize the above local stability of all equilibria in the following proposition.

Proposition 2.1.

(i)

The trivial equilibrium $E_{0}$ is unstable.
(ii)

The semi-trivial equilibrium $E_{1}$ is unstable.
(iii)

The semi-trivial equilibrium $E_{2}$ is globally asymptotically stable if (H1) does not hold.
(iv)

The semi-trivial equilibrium $E_{12}$ is stable if $\mu>a_{31}u^{*}_{12}$ .
(v)

The semi-trivial equilibrium $E_{13}$ is unstable.
(vi)

The positive equilibrium $E_{*}$ exists and is stable if (H3) holds.

Furthermore, we can obtain the following two global results.

Theorem 2.1.

Let assumption (H1) and $\mu\geq a_{31}u^{*}_{12}$ hold. Then the positive equilibria $E_{12}$ exists, and it is globally asymptotically stable.

Proof.

Modify the standard Lyapunov function in the following form

V(u(t),v(t),w(t))=\frac{a_{21}}{a_{12}}(u-u^{*}_{12}\ln u)+v-v^{*}_{12}\ln v+\frac{a_{13}a_{21}}{a_{12}a_{31}}w.

Then, with equation (2.2),

	$\displaystyle\frac{d}{dt}$	$\displaystyle V(u(t),v(t),w(t))$
		$\displaystyle=\frac{a_{21}}{a_{12}}(u-u^{*}_{12})(r_{1}(1-u)-a_{12}v-a_{13}w)+$
		$\displaystyle\quad(v-v^{*}_{12})(r_{2}(1-v)+a_{21}u)+\frac{a_{13}a_{21}}{a_{12}a_{31}}(-\mu w+a_{31}u)$
		$\displaystyle=-\frac{r_{1}a_{21}}{a_{12}}(u-u^{}_{12})^{2}-r_{2}(v-v^{}_{12})^{2}+\frac{a_{13}a_{21}}{a_{12}}(\frac{-\mu}{a_{31}}+u^{*}_{12})w\leq 0.$

Hence $\{dV/dt=0\}=\{(u^{*}_{12},v^{*}_{12},0)\}$ if $\mu>a_{31}u^{*}_{12}$ , or $\{dV/dt=0\}=\{W\geq 0:(u^{*}_{12},v^{*}_{12},W)\}$ if $\mu=a_{31}u^{*}_{12}$ . However, the maximal invariant set of $\{dV/dt=0\}$ is the singleton set $\{(u^{*}_{12},v^{*}_{12},0)\}$ for these two possibilities. By LaSalle’s Invariance Principle, we show the global stability of equilibrium $E_{12}$ if $\mu\geq a_{31}u^{*}_{12}$ . The proof is completed.

∎

Theorem 2.2.

Let assumption (H3) hold. Then the positive equilibria $E_{*}$ exists, and it is globally asymptotically stable.

Proof.

Define the Lyapunov function in the following form

V(u(t),v(t),w(t))=\frac{a_{12}}{a_{21}}(u-u_{*}\ln u)+v-v_{*}\ln v+\frac{a_{13}a_{21}}{a_{12}a_{31}}(w-w_{*}\ln w).

Then

	$\displaystyle\frac{d}{dt}$	$\displaystyle V(u(t),v(t),w(t))$
	$\displaystyle=$	$\displaystyle\frac{a_{12}}{a_{21}}(u-u_{})(r_{1}(1-u)-a_{12}v-a_{13}w)+(v-v_{})(r_{2}(1-v)+a_{21}u)$
		$\displaystyle+\frac{a_{13}a_{21}}{a_{12}a_{31}}(w-w_{*})(-\mu+a_{31}u)$
	$\displaystyle=$	$\displaystyle\frac{a_{12}}{a_{21}}(u-u_{})(r_{1}(u_{}-u)+a_{12}(v_{}-v)+a_{13}(w_{}-w))$
		$\displaystyle+(v-v_{})(r_{2}(v_{}-v)+a_{21}(u-u_{}))+\frac{a_{13}a_{21}}{a_{12}a_{31}}(w-w_{})(a_{31}(u-u_{*}))$
	$\displaystyle=$	$\displaystyle-\frac{r_{1}a_{21}}{a_{12}}(u-u_{})^{2}-r_{2}(v-v_{})^{2}\leq 0.$

Hence by the LaSalle’s Invariance Principle, the $\omega$ -limit set of any solution of (1.2) is contained in the maximal invariant subset of $\{dV/dt=0\}=\{(u_{*},v_{*},W):W>0\}$ , which is the singleton $\{(u_{*},v_{*},w_{*})\}$ . We complete the proof. ∎

3 Existence of traveling wave solutions

In this section, motivated by Huang [10, 11], the high dimensional shooting method is implemented to investigate the existence of wave fronts, or traveling wave solutions of (1.1) from $E_{1}$ to $E_{*}$ . Following the ideas of Huang, we list the main steps as follows.

(i)

By using the moving coordinates, the reaction-diffusion system is transformed into a ODE system.
(ii)

Construct a variant of Wazewski set $\Sigma$ with $E_{1}$ as its boundary point and also containing $E_{*}$ . Dynamics of the system on all boundaries of $\Sigma$ should be clarified.
(iii)

Analyze the structure of unstable manifold of $E_{1}$ , and pick up a curve contained in the unstable manifold of $E_{1}$ with two end points on the “exit set” of boundary of $\Sigma$ .
(iv)

Show that there exists a particular point on this curve, and the solution starting from this point will stay in the interior of $\Sigma$ for all $t\geq 0$ .
(v)

Define a nonempty subset of $\Sigma$ which contains the point stayed in $\Sigma$ for all positive time under the action of the ODE system. Then by constructing a Lyapunov function and using the LaSalle’s invariance principle, we obtain the main result.

3.1 The ODE forms and the Lienard Transformation

We consider the solution of (1.1) with the moving coordinate $\xi$ and wave speed $c$ of the form

\left\{\begin{aligned} u(x,t)=U(x+ct)=U(\xi),\\ v(x,t)=V(x+ct)=V(\xi),\\ w(x,t)=W(x+ct)=W(\xi),\end{aligned}\right.

(3.1)

satisfying the asymptotical boundary conditions from $E_{1}$ to $E_{*}$ , that is,

	$\displaystyle\lim_{\xi\to-\infty}\big{(}U(\xi),V(\xi),W(\xi)\big{)}$	$\displaystyle=(1,0,0),\quad\text{ and }$		(3.2)
	$\displaystyle\lim_{\xi\to\infty}\big{(}U(\xi),V(\xi),W(\xi)\big{)}$	$\displaystyle=(u_{},v_{},w_{*}).$		(3.2)

A direct computation shows that $(U(\xi),V(\xi),W(\xi))$ is a traveling wave solution of (1.1) if and only if $(U(\xi),V(\xi),W(\xi))$ is a solution of the system,

\left\{\begin{aligned} cU^{\prime}&=U\left(r_{1}(1-U)-a_{12}V-a_{13}W\right),\\ cV^{\prime}&=V(r_{2}(1-V)+a_{21}U),\\ cW^{\prime}&=dW^{\prime\prime}+W(-\mu+a_{31}U).\end{aligned}\right.

(3.3)

For a constant $c>0$ we consider the Lienard transformation to make the following changes of variables and scaling:

	$\displaystyle X_{1}(t)$	$\displaystyle=U(ct),$
	$\displaystyle X_{2}(t)$	$\displaystyle=V(ct),$
	$\displaystyle Y(t)$	$\displaystyle=W(ct),$
	$\displaystyle Z(t)$	$\displaystyle=\frac{1}{c}[cW(ct)-dW^{\prime}(ct)].$

Then, upon straightforward computation, (3.3) is transformed to a four dimensional system,

$\displaystyle\dot{X}_{1}$	$\displaystyle=X_{1}\left(r_{1}(1-X_{1})-a_{12}X_{2}-a_{13}Y\right),$	(3.4)
$\displaystyle\dot{X}_{2}$	$\displaystyle=X_{2}(r_{2}(1-X_{2})+a_{21}X_{1}),$
$\displaystyle\dot{Y}$	$\displaystyle=\rho[Y-Z],\ \ \ \ \rho=\frac{c^{2}}{d},$
$\displaystyle\dot{Z}$	$\displaystyle=Y(-\mu+a_{31}X_{1}).$

It is clear that $(U(\xi),V(\xi),W(\xi))$ is a nonnegative solution of (3.3) with asymptotical boundary conditions (3.2) if and only if $\left(X_{1}(t),X_{2}(t),Y(t),Z(t)\right)$ is a nonnegative solution of (3.4) satisfying the asymptotical boundary conditions,

\left(X_{1}(-\infty),X_{2}(-\infty),Y(-\infty),Z(-\infty)\right)=(1,0,0,0),

\left(X_{1}(\infty),X_{2}(\infty),Y(\infty),Z(\infty)\right)=(u_{*},v_{*},w_{*},w_{*}).

3.2 The Wazewski set and its Exist Subset

Let $\sigma_{1}$ and $\sigma_{2}$ be constants defined by

\sigma_{1}=\frac{\rho+\sqrt{\rho^{2}-4\rho(a_{31}-\mu)}}{2\rho}=\frac{c+\sqrt{c^{2}-4d(a_{31}-\mu)}}{2c}

and

\sigma_{2}=\frac{\rho+\sqrt{\rho^{2}+4\rho\mu}}{2\rho}.

It is clear that $0<\sigma_{1}<1$ and $\sigma_{2}>1$ if $c>c_{*}\equiv 2\sqrt{d(a_{31}-\mu)}$ or, equivalently, $\rho>4(a_{31}-\mu)$ . For $c>c_{*}$ , we define a wedged like region $\Sigma\subset\mathbb{R}^{4}$ as follows,

\Sigma=\{(X_{1},X_{2},Y,Z):0\leq X_{1}\leq 1,0\leq X_{2}\leq 1+\frac{a_{21}}{r_{2}},Y\geq 0,\sigma_{1}Y\leq Z\leq\sigma_{2}Y\}.

Then the boundary of $\Sigma$ consists of surfaces $P_{1}$ , $P_{2}$ , $Q_{1}\sim Q_{5}$ represented by

	$\displaystyle Q_{1}$	$\displaystyle=\{X_{1}=0,0\leq X_{2}\leq 1+\frac{a_{21}}{r_{2}},Y\geq 0,\sigma_{1}Y\leq Z\leq\sigma_{2}Y\},$
	$\displaystyle Q_{2}$	$\displaystyle=\{X_{1}=1,0\leq X_{2}\leq 1+\frac{a_{21}}{r_{2}},Y\geq 0,\sigma_{1}Y\leq Z\leq\sigma_{2}Y\},$
	$\displaystyle Q_{3}$	$\displaystyle=\{0\leq X_{1}\leq 1,X_{2}=0,Y\geq 0,\sigma_{1}Y\leq Z\leq\sigma_{2}Y\},$
	$\displaystyle Q_{4}$	$\displaystyle=\{0\leq X_{1}\leq 1,X_{2}=1+\frac{a_{21}}{r_{2}},Y\geq 0,\sigma_{1}Y\leq Z\leq\sigma_{2}Y\},$
	$\displaystyle Q_{5}$	$\displaystyle=\{0<X_{1}<1,0<X_{2}<1+\frac{a_{21}}{r_{2}},Y=Z=0\},$
	$\displaystyle P_{1}$	$\displaystyle=\{0<X_{1}<1,0<X_{2}<1+\frac{a_{21}}{r_{2}},Y>0,Z=\sigma_{1}Y\},$
	$\displaystyle P_{2}$	$\displaystyle=\{0<X_{1}<1,0<X_{2}<1+\frac{a_{21}}{r_{2}},Y>0,Z=\sigma_{2}Y\}.$

The vector field of (3.4) has a very simple property in the boundary of $\Sigma$ , which can be characterized by the following two lemmas.

Lemma 3.1.

Let $c>c_{*}$ and $\Phi_{t}(p)$ be the flow of (3.4), i.e. $\Phi_{t}(p)$ is a solution of (3.4) satisfying the initial condition $\Phi_{0}(p)=p\in\mathbb{R}^{4}$ . Then for any $p\in Int(\Sigma)$ , $\Phi_{t}(p)$ cannot leave $\Sigma$ from a point in the boundary $\cup_{i=1}^{5}Q_{i}$ of $\Sigma$ at any positive time.

Proof.

It is obvious that the boundary $Q_{1}$ , $Q_{3}$ and $Q_{5}$ are invariant sets of (3.4). Hence any solution of (3.4) through a point in the interior of $\Sigma$ can not leave $\Sigma$ from a point in the set $Q_{1}\cup Q_{3}\cup Q_{5}$ .

First suppose $p_{0}=(1,X_{2}(0),Y(0),Z(0))\in Q_{2}$ , then the first equation of (3.4) yields that at $p_{0}$ ,

\dot{X}_{1}(0)=-a_{12}X_{2}(0)-a_{13}Y(0)<0.

So that the vector field of (3.4) points interior of $\Sigma$ and hence $\Phi_{t}(p_{0})$ cannot exit $\Sigma$ at $p_{0}$ in the set $Q_{2}$ .

Next suppose $p_{1}=(X_{1}(0),1+\frac{a_{21}}{r_{2}},Y(0),Z(0))\in Q_{4}$ , then the second equation of (3.4) yields that at $p_{1}$ ,

	$\displaystyle\dot{X}_{2}(0)$	$\displaystyle=\left(1+\frac{a_{21}}{r_{2}}\right)\left[r_{2}\left(1-\left(1+\frac{a_{21}}{r_{2}}\right)\right)+a_{21}X_{1}(0)\right]$
		$\displaystyle=-a_{21}-\frac{a_{21}^{2}}{r_{2}}+\left(a_{21}+\frac{a_{21}^{2}}{r_{2}}\right)X_{1}(0)$
		$\displaystyle<-a_{21}-\frac{a_{21}^{2}}{r_{2}}+a_{21}+\frac{a_{21}^{2}}{r_{2}}=0.$

So that the vector field of (3.4) points interior of $\Sigma$ and hence $\Phi_{t}(p_{1})$ cannot exit $\Sigma$ at $p_{1}$ in the set $Q_{4}$ . ∎

Next let us study the vector field at the sets $P_{1}$ and $P_{2}$ .

Lemma 3.2.

Let $c>c_{*}$ , where $c$ is the wave speed in the system (3.4). Then the vector field of (3.4) at $p\in P_{1}\cup P_{2}$ points outside of $\Sigma$ .

Proof.

First consider a point $p_{1}=(X_{1}(0),X_{2}(0),Y(0),Z(0))\in P_{1}$ . Then $Z(0)=\sigma_{1}Y(0)<Y(0)$ and $Y(0)>0$ . By the last two equations of (3.4), we obtain

\dot{Y}(0)=\rho[Y(0)-Z(0)]>0,

	$\displaystyle\frac{\dot{Z}(0)}{\dot{Y}(0)}$	$\displaystyle=\frac{Y(0)(-\mu+a_{31}X_{1}(0))}{\rho[Y(0)-Z(0)]}$
		$\displaystyle=\frac{-\mu+a_{31}X_{1}(0)}{\rho(1-\sigma_{1})}$
		$\displaystyle<\frac{a_{31}-\mu}{\rho(1-\sigma_{1})}=\sigma_{1}.$

It implies that the vector field at the point $p_{1}\in P_{1}$ points outside of $\Sigma$ . Next, let $p_{2}=(X_{1}(0),X_{2}(0),Y(0),Z(0))\in P_{2}$ . Then $Z(0)=\sigma_{2}Y(0)>Y(0)$ and $Y(0)>0$ . At $p_{2}$ , we have

\dot{Y}(0)=\rho[Y(0)-Z(0)]<0,

	$\displaystyle\frac{\dot{Z}(0)}{\dot{Y}(0)}$	$\displaystyle=\frac{Y(0)(-\mu+a_{31}X_{1}(0))}{\rho[Y(0)-Z(0)]}$
		$\displaystyle=\frac{-a_{31}X_{1}(0)+\mu}{\rho(\sigma_{2}-1)}$
		$\displaystyle<\frac{\mu}{\rho(\sigma_{2}-1)}=\sigma_{2}.$

It implies that the vector field at the point $p_{2}\in P_{2}$ points outside of $\Sigma$ . ∎

By Lemma 3.1 and Lemma 3.2, we can see that with initial point $p\in\Sigma$ $\Phi_{t}(p)$ can only leave $\Sigma$ from a point in $P_{1}\cup P_{2}$ at some positive time.

3.3 The Unstable Manifold of $E_{1}$

Now we turn to study the unstable manifold of the equilibrium $E_{1}$ . The linearized system of (3.4) at $E_{1}$ is

$\displaystyle\dot{X_{1}}$	$\displaystyle=-r_{1}X_{1}-a_{12}X_{2}-a_{13}Y,$	(3.5)
$\displaystyle\dot{X_{2}}$	$\displaystyle=(r_{2}+a_{21})X_{2},$
$\displaystyle\dot{Y}$	$\displaystyle=\rho[Y-Z]$
$\displaystyle\dot{Z}$	$\displaystyle=(-\mu+a_{31})Y,$

with Jacobian matrix of (3.4) at $E_{1}$ ,

\begin{bmatrix}-r_{1}&-a_{12}&-a_{13}&0\\ 0&r_{2}+a_{21}&0&0\\ 0&0&\rho&-\rho\\ 0&0&-\mu+a_{31}&0\end{bmatrix}.

For $c>c_{*}$ , upon a direct computation, the matrix has one negative eigenvalue

\lambda_{0}=-r_{1}

and three distinct positive eigenvalues

$\displaystyle\lambda_{1}$	$\displaystyle=r_{2}+a_{21},$	(3.6)
$\displaystyle\lambda_{2}$	$\displaystyle=\frac{\rho+\sqrt{\rho^{2}-4\rho(a_{31}-\mu)}}{2},$
$\displaystyle\lambda_{3}$	$\displaystyle=\frac{\rho-\sqrt{\rho^{2}-4\rho(a_{31}-\mu)}}{2},$

where the corresponding eigenvectors to $\lambda_{1}$ , $\lambda_{2}$ and $\lambda_{3}$ are

$\displaystyle h_{1}$	$\displaystyle=\left[-\frac{a_{12}}{\lambda_{1}+r_{1}},1,0,0\right]^{T},$	(3.7)
$\displaystyle h_{2}$	$\displaystyle=\left[-\frac{a_{13}}{\lambda_{2}+r_{1}},0,1,\frac{a_{31}-\mu}{\lambda_{2}}\right]^{T},$
$\displaystyle h_{3}$	$\displaystyle=\left[-\frac{a_{13}}{\lambda_{3}+r_{1}},0,1,\frac{a_{31}-\mu}{\lambda_{3}}\right]^{T}.$

Hence, the unstable manifold $\varepsilon_{1}^{u}$ of $E_{1}$ is tangent to the plane

P=\{k_{1}h_{1}+k_{2}h_{2}+k_{3}h_{3}+E_{1}:k_{1},k_{2},k_{3}\in\mathbb{R}\}=\{x_{1},x_{2},y,z\}

where

	$\displaystyle x_{1}$	$\displaystyle=1-\frac{k_{1}a_{12}}{\lambda_{1}+r_{1}}-\frac{k_{2}a_{13}}{\lambda_{2}+r_{1}}-\frac{k_{3}a_{13}}{\lambda_{3}+r_{1}},$
	$\displaystyle x_{2}$	$\displaystyle=k_{1},$
	$\displaystyle y$	$\displaystyle=k_{2}+k_{3},$
	$\displaystyle z$	$\displaystyle=\frac{k_{2}(a_{31}-\mu)}{\lambda_{2}}+\frac{k_{3}(a_{31}-\mu)}{\lambda_{3}}.$

Lemma 3.3.

There are two points $p_{1}^{*}\in P_{1}$ and $p_{2}^{*}\in P_{2}$ and a curve $\gamma\in\varepsilon_{1}^{u}$ which connects $p_{1}^{*}$ and $p_{2}^{*}$ such that

\gamma\setminus\{p_{1}^{*},p_{2}^{*}\}\subset{\rm Int}(\Sigma).

Proof.

Note that the transformations: $(x_{2},y,z)\mapsto(k_{1},k_{2},k_{3})$ ,

	$\displaystyle x_{2}$	$\displaystyle=k_{1},$
	$\displaystyle y$	$\displaystyle=k_{2}+k_{3},$
	$\displaystyle z$	$\displaystyle=\frac{k_{2}(a_{31}-\mu)}{\lambda_{2}}+\frac{k_{3}(a_{31}-\mu)}{\lambda_{3}}$

is invertible. Thus the plane can also be expressed as

P=\left\{\left(1-\frac{x_{2}a_{12}}{\lambda+r_{2}}-\frac{\lambda_{2}\lambda_{3}z+r_{1}(a_{31}-\mu)y}{(a_{31}-\mu)(\lambda_{2}+r_{1})(\lambda_{3}+r_{1})},x_{2},y,z\right):x_{2},y,z\in\mathbb{R}\right\}.

F Since the unstable manifold $\varepsilon_{1}^{u}$ is tangent to $P$ at $E_{1}$ , there is a small $k>0$ such that the unstable manifold $\varepsilon_{1}^{u}$ of $E_{1}$ can be expressed as

\varepsilon_{1}^{u}=\left\{\left(x_{1}(x_{2},y,z),x_{2},y,z\right):(x_{2},y,z)\in[0,k]^{3}\right\}

where

\displaystyle x_{1}(x_{2},y,z)=1-\frac{x_{2}a_{12}}{\lambda+r_{2}}-\frac{\lambda_{2}\lambda_{3}z+r_{1}(a_{31}-\mu)y}{(a_{31}-\mu)(\lambda_{2}+r_{1})(\lambda_{3}+r_{1})}+\eta(x_{2},y,z),

(3.8)

and the function $\eta(x_{2},y,z)$ satisfies

\displaystyle\eta(x_{2},y,z)=O(x_{2}^{2}+y^{2}+z^{2})\ \ as\ \ (x_{2},y,z)\rightarrow(0,0,0).

(3.9)

Select a sufficiently small $\epsilon>0$ with $k\geq\sigma_{2}\epsilon$ . It follows from (3.8) and (3.9) that

\displaystyle 0<x_{1}(x_{2},y,z)<1,\ \ for\ x_{2}\in[0,\epsilon],\ y\in[0,\epsilon],\ z\in[\sigma_{1}y,\sigma_{2}y].

(3.10)

Now, we define $\gamma\subset\varepsilon_{1}^{u}$ as

\gamma=\left\{(x_{1}(\epsilon,\epsilon,z),\epsilon,\epsilon,z):z\in[\sigma_{1}\epsilon,\sigma_{2}\epsilon]\right\},

and let $p_{1}^{*}=(x_{1}(\epsilon,\epsilon,\sigma_{1}\epsilon),\epsilon,\epsilon,\sigma_{1}\epsilon)$ and $p_{2}^{*}=(x_{1}(\epsilon,\epsilon,\sigma_{2}\epsilon),\epsilon,\epsilon,\sigma_{2}\epsilon)$ . Then it is clear that $p_{i}^{*}\in P_{i}$ , for $i=1,2$ and $\gamma$ is a curve connecting $p_{1}^{*}$ and $p_{2}^{*}$ . Moreover, by (3.10), it is obvious that $\gamma\backslash\{p_{1}^{*},p_{2}^{*}\}\subset{\rm Int}({\Sigma})$ which is denoted the interior set of $\Sigma$ . ∎

Lemma 3.4.

For each $c\geq c_{*}$ , there is a point $p_{*}\in\varepsilon_{1}^{u}\cap{\rm Int}(\Sigma)$ such that the solution $\Phi_{t}(p_{*})$ of (3.4) through the $p_{*}$ stays in ${\rm Int}(\Sigma)$ for all $t\geq 0$ .

Proof.

Let $\gamma\subset\varepsilon_{1}^{u}\cap{\rm Int}(\Sigma)$ be defined as in Lemma 3.3. We define subsets $\gamma_{1}$ and $\gamma_{2}$ of $\gamma$ as follows:

\gamma_{i}=\{p\in\gamma:\ there\ is\ a\ t\geq 0\ such\ that\ \Phi_{t}(p)\in P_{i}\}\ \ for\ i=1,\ 2.

By Lemma 3.3, it is clear that $\gamma_{1}$ and $\gamma_{2}$ are nonempty, and by the Lemma 3.2, the vector field of (3.4) at each point in the plane $P_{i}(i=1,\ 2)$ points to the exterior of $\Sigma$ if $c\geq c_{*}$ . By the continuity of solutions on initial values, we deduce that two sets $\gamma_{1}$ and $\gamma_{2}$ are open relative to $\gamma$ . Hence, $\gamma\backslash(\gamma_{1}\cup\gamma_{2})\neq\emptyset$ , because of connectedness of $\gamma$ . Let $p_{*}\in\gamma\backslash(\gamma_{1}\cup\gamma_{2})$ . Then the definitions of $\gamma$ and $\gamma_{i}$ imply that

\displaystyle\Phi_{t}(p_{*})\in{\rm Int}(\Sigma)\ \ for\ all\ \ t\geq 0.

∎

3.4 Existence and Non-existence of Traveling wave solutions from $E_{1}$ to $E_{*}$

Let us define an essential subsets of $\Sigma$ :

G=\{p\in{\rm Int}(\Sigma):\Phi_{p}(t)\in{\rm Int}(\Sigma)\ for\ all\ t\geq 0\}.

By Lemma 3.4, we know that $G$ is nonempty. The desired heteroclitic orbit from $E_{1}$ to $E_{*}$ will be obtained in $G$ .

Proposition 3.1.

Let assumption $(H_{3})$ hold. Then system (1.1) has a traveling wave solution connecting $E_{1}$ and $E_{*}$ for wave speed $c>c_{*}$ .

Proof.

Define the Lyapunov function $L:[0,\infty)\rightarrow\mathbb{R}$ , which is associated solutions of (3.4) with initial in $G$ ,

	$\displaystyle L(t)=$	$\displaystyle\frac{a_{21}}{a_{12}}\left(X_{1}-u_{}\ln X_{1}\right)+\left(X_{2}-v_{}\ln X_{2}\right)$
		$\displaystyle+\frac{a_{13}a_{21}}{a_{12}a_{31}}\left(Y-w_{*}\ln Y-(Y-Z)(1-\frac{1}{Y})\right),$

and it is clear that $L$ is well defined and continuously differentiable on $G$ . Moreover, the orbital derivative of $L(t)$ along (3.4) is

$\displaystyle\dot{L}(t)=$	$\displaystyle\frac{a_{21}}{a_{12}}\left(\dot{X_{1}}-u_{}\frac{\dot{X_{1}}}{X_{1}}\right)+\left(\dot{X_{2}}-v_{}\frac{\dot{X_{2}}}{X_{2}}\right)+\frac{a_{13}a_{21}}{a_{12}a_{31}}\left(\dot{Z}-w_{*}\frac{\dot{Z}}{Z}-\rho\frac{(Y-Z)^{2}}{Y^{2}}\right)$	(3.11)
$\displaystyle=$	$\displaystyle\frac{a_{21}}{a_{12}}(X_{1}-u_{})\left(r_{1}(u_{}-X_{1})+a_{12}(v_{}-X_{2})+a_{13}(w_{}-Y)\right)$
	$\displaystyle+(X_{2}-v_{})\left(r_{2}(v_{}-X_{2})+a_{21}(X_{1}-u_{*})\right)$
	$\displaystyle+\frac{a_{13}a_{21}}{a_{12}a_{31}}(Y-w_{})\left(a_{13}(X_{1}-u_{})\right)-\frac{\rho a_{13}a_{21}}{a_{12}a_{31}}\frac{(Y-Z)^{2}}{Y^{2}}$
$\displaystyle=$	$\displaystyle-\frac{a_{21}r_{1}}{a_{12}}(X_{1}-u_{})^{2}-r_{2}(X_{2}-v_{})^{2}-\frac{\rho a_{13}a_{21}}{a_{12}a_{31}}\frac{(Y-Z)^{2}}{Y^{2}}\leq 0.$

Hence, by the LaSalle’s Invariance Principle, the $\omega$ -limit set of any solution of (3.4) is contained in the largest invariant subset of $\{dL/dt=0\}=\{(u_{*},v_{*},W,W):W>0\}$ , which is the singleton $\{(u_{*},v_{*},w_{*},w_{*})\}$ . This completes the proof. ∎

Proposition 3.2.

The system (1.1) does not have a positive traveling wave solution connecting $E_{1}$ and $E_{*}$ for wave speed $0<c<c_{*}$ .

Proof.

If $0<c<c_{*}$ , then $\lambda_{2}$ and $\lambda_{3}$ in (3.6) are a pair of complex eigenvalues with positive real parts. By looking at the associated eigenvectors given in (3.7), one is able to conclude that if a solution is in the unstable manifold of the equilibrium $E_{1}$ , then its $w(t)$ component can not keep nonnegative all the time when the solution converge to $E_{1}$ as $t\to-\infty$ . This completes the proof. ∎

4 Numerical Simulations and Brief Discussions

We close this work by performing some interesting numerical simulations with FiPy which is an object oriented PDE solver, written in Python [5]. Some brief discussions and biological implications are also given.

4.1 Numerical Simulations

Consider system (1.1) in one-dimensional bounded space $[0,10]$ with Neumann boundary condition and parameters, $r_{1}=0.7$ , $r_{2}=0.3$ , $\mu=0.15$ , $a_{12}=0.15$ , $a_{13}=0.5$ , $a_{21}=0.2$ , $a_{31}=0.5$ . It is easy to check that the assumptions (H1)-(H3) are satisfied for these parameters. Three cases of initial functions are taken in the following.

(i)

Let $u(x,0)=0.2$ and $v(x,0)=0.1$ for all $x$ , and

$w(x,0)=\left\{\begin{aligned} 0.08&&\text{ if }4.8\leq x\leq 5.2,\\ 0&&\text{ otherwise.}\end{aligned}\right.$
(ii)

Let $u(x,0)=1$ for all $x$ ,

$v(x,t)=\left\{\begin{aligned} 0.1&&\text{ if }0\leq x\leq 5,\\ 0&&\text{ otherswise.}\end{aligned}\right.$

$w(x,0)=\left\{\begin{aligned} 0.1&&\text{ if }5\leq x\leq 10,\\ 0&&\text{ otherwise.}\end{aligned}\right.$
(iii)

Let $v(x,0)=1$ for all $x$ ,

$u(x,t)=\left\{\begin{aligned} 0.15&&\text{ if }0\leq x\leq 5,\\ 0&&\text{ otherwise.}\end{aligned}\right.$

$w(x,0)=\left\{\begin{aligned} 0.1&&\text{ if }0\leq x\leq 5,\\ 0&&\text{ otherwise.}\end{aligned}\right.$

Please refer Figure 1-3(a) for graphical representations of initial functions.

For the first case, we simulate that three species $u$ , $v$ , and $w$ migrate into a new habitat simultaneously. The initial condition is close to $E_{0}$ at the time $t=0$ . Note that species $w$ appear in some places of the habitat for $t=0$ , and it diffuses to the whole habitat as time increasing. Meanwhile, species $u$ is affected by the diffusion of species $w$ negatively, and species $v$ is affected positively by $u$ . So, eventually, all population densities of all species approach to the positive spatial homogeneous state solution $E_{*}$ of (1.1). We can see easily the process of transition from $E_{0}$ to $E_{*}$ presented in Figure 1 (b).

For case 2, different initial functions which are close to $E_{1}$ at time $t=0$ are considered. We would like to simulate that species $v$ and $w$ invade simultaneously into a new habitat with saturated species $u$ . Note that species $v$ cannot diffuse, so species $v$ appears in the left half space showed in Figure 2(b) for all time. On the left half space, the transition process from $E_{1}$ to $E_{*}$ can be observed clearly. However, on the right half space, we conjecture that the transition process from $E_{1}$ to $E_{13}$ happens.

For the last case, we set the initial functions close to $E_{2}$ which is saturated by species $v$ . And species $u$ and $w$ migrate to the left half space of the habitat with little populations at time $t=0$ . The transition process from $E_{2}$ to $E_{*}$ on the left half space can also be observed. On the right half space of habitat, $u=w=0$ for time $t=0$ . However, species $w$ can diffuse to the right half space of habitat. We conjecture that the population of $w$ in the right half space of habitat is proportional to the diffusion coefficient $d$ .

4.2 Brief Discussions

In this work, for system without diffusion (1.2), we have showed that if (H1) doest not hold then $\lim_{t\to\infty}u(t)=0$ and $\lim_{t\to\infty}w(t)=0$ (Lemma 2.2). It can be showed that system (1.2) asymptotically approaches to a one-dimensional subsystem only involving species $v$ , and this implies that $E_{2}$ is globally asymptotically stable by Markus limiting theorem [16]. If (H2) doest not hold then we only have $\lim_{t\to\infty}w(t)=0$ (Lemma 2.3). However, we obtain a global results of $E_{12}$ in Theorem 2.1 under the assumptions (H1) and $\mu\geq a_{31}u^{*}_{12}$ . It is clearly that $\mu\geq a_{31}$ is a sufficient condition of $\mu\geq a_{31}u^{*}_{12}$ because of $0<u^{*}_{12}<1$ . Finally, the assumption (H3) is a sufficient and necessary condition for the existence and global stability of the positive equilibrium $E_{*}$ (Thoerem 2.2). Since the positive equilibrium $E_{*}$ does not exist if (H3) is not true, and condition (H3) can be rewritten as the form,

0<r_{1}r_{2}a_{31}-r_{1}r_{2}\mu-a_{12}a_{31}r_{2}-a_{12}a_{21}\mu=(r_{1}r_{2}+a_{12}a_{21})(u^{*}_{12}a_{31}-\mu).

All global dynamics with respective to essential parameters are summarized in Table 1.

	$E_{2}$	$E_{12}$	$E_{*}$
$r_{1}\leq a_{12}$	GAS	$\nexists$	$\nexists$
$r_{1}>a_{12}$ and $\mu\geq a_{31}u^{*}_{12}$	Saddle	GAS	$\nexists$
$r_{1}>a_{12}$ and $\mu<a_{31}u^{*}_{12}$	Saddle	Saddle	GAS

Table 1: Classification of dynamics of equilibria with respective to all essential parameters

For system with diffusion (1.1) and using the high-dimensional shooting method, we have show the existence of transition between two different spatial homogeneous states, that is, the existence of traveling wave solution from $E_{1}$ to $E_{*}$ for $c>c_{*}=2\sqrt{d(a_{31}-\mu)}$ which is called the minimal speed. By the similar manner, motivated by Table 1, it is naturally to consider the possibilities of other cases. For example, the traveling wave solution from $E_{2}$ or $E_{12}$ to $E_{*}$ under the assumptions $r_{1}>a_{12}$ and $\mu<a_{31}u^{*}_{12}$ . Or the traveling wave solution from $E_{1}$ or $E_{2}$ to $E_{12}$ under the assumptions $r_{1}>a_{12}$ and $\mu\geq a_{31}u^{*}_{12}$ . We left these open problems to be considered.

References

[1] X. Chen Existence, uniqueness, and asymptotic stability of traveling waves in nonlocal evolution equations, Advances in Differential Equations, 2(1) : 125–160, 1997.
[2] S. R. Dunbar. Travelling wave solutions of diffusive Lotka-Volterra equations. Journal of Mathematical Biology, 17(1):11–32, 1983.
[3] S. R. Dunbar. Traveling wave solutions of diffusive Lotka-Volterra equations: a heteroclinic connection in $\mathbb{R}^{4}$ . Transactions of the American Mathematical Society, 557–594, 1984.
[4] S. R. Dunbar. Traveling waves in diffusive predator-prey equations: periodic orbits and point-to-periodic heteroclinic orbits. SIAM Journal on Applied Mathematics, 46(6):1057–1078, 1986.
[5] J. E. Guyer, D. Wheeler and J. A. Warren. FiPy: Partial Differential Equations with Python Computing in Science & Engineering, 11(3):6–15, 2009.
[6] R. A. Fisher. The wave of advance of advantageous genes. Annals of eugenics, 7(4):355–369, 1937.
[7] N. J. Ford, P. M. Lumb, and E. Ekaka-a. Mathematical modelling of plant species interactions in a harsh climate. Journal of computational and Applied Mathematics, 234(9):2732–2744, 2010.
[8] R. Gardner and J. A. Smoller. The existence of periodic travelling waves for singularly perturbed predator-prey equations via the conley index. 1983.
[9] C. Hsu, C. Yang, T. Yang, and T. Yang. Existence of traveling wave solutions for diffusive predator-prey type systems. Journal of Differential Equations, 252(4):3040–3075, 2012.
[10] W. Huang. Traveling wave solutions for a class of predator-prey systems. Journal of Dynamics and Differential Equations, 24(3):633–644, 2012.
[11] W. Huang. A geometric approach in the study of traveling waves for some classes of non-monotone reaction-diffusion systems. Journal of Differential Equations, 260(3):2190–2224, 2016.
[12] C. K. Jones, R. Gardner, and T. Kapitula. Stability of travelling waves for non-convex scalar viscous conservation laws. Communications on Pure and Applied Mathematics, 46(4):505–526, 1993.
[13] P. I. Kolmogorov, A. and N. Piskunov. A study of the equation of diffusion with increase in the quantity of matter, and its application to a biological problem. Moscow University Mathematics Bulletin, 1:1–26, 1937.
[14] Liang X, Zhao XQ. Asymptotic speeds of spread and traveling waves for monotone semiflows with applications. Communications on pure and applied mathematics 2007;60:1-40.
[15] J. Lin, W. Wang, C. Zhao, and T. Yang. Global dynamics and traveling wave solutions of two predators-one prey models. Discrete & Continuous Dynamical Systems-Series B, 20(4), 2015.
[16] L. Markus. Asymptotically autonomous differential systems,. Contributions to the Theory of Nonlinear Oscillations, 3:17–29, 1956.
[17] R. McGehee and R. A. Armstrong. Some mathematical problems concerning the ecological principle of competitive exclusion. Journal of Differential Equations, 23(1):30–52, 1977.
[18] J. Murray. Mathematical Biology II Spatial Models and Biomedical Applications. Springer, 2003.
[19] T. Namba, K. Tanabe, and N. Maeda. Omnivory and stability of food webs. Ecological Complexity, 5(2):73–85, 2008.
[20] H. L Smith and X.-Q. Zhao. Global asymptotic stability of traveling waves in delayed reaction-diffusion equations. SIAM Journal on Mathematical Analysis, 31(3) : 514–534, 2000.
[21] A. I. Volpert, V. A. Volpert, and V. A. Volpert. Traveling wave solutions of parabolic systems. American Mathematical Society, 1994.
[22] X. Zou and J. Wu. Existence of traveling wave fronts in delayed reaction-diffusion systems via the monotone iteration method Proceedings of the American Mathematical Society, 125(9):2589–2598, 1997.

$\displaystyle\dot{L}(t)=$	$\displaystyle\frac{a_{21}}{a_{12}}\left(\dot{X_{1}}-u_{}\frac{\dot{X_{1}}}{X_{1}}\right)+\left(\dot{X_{2}}-v_{}\frac{\dot{X_{2}}}{X_{2}}\right)+\frac{a_{13}a_{21}}{a_{12}a_{31}}\left(\dot{Z}-w_{*}\frac{\dot{Z}}{Z}-\rho\frac{(Y-Z)^{2}}{Y^{2}}\right)$	(3.11)
$\displaystyle=$	$\displaystyle\frac{a_{21}}{a_{12}}(X_{1}-u_{})\left(r_{1}(u_{}-X_{1})+a_{12}(v_{}-X_{2})+a_{13}(w_{}-Y)\right)$
	$\displaystyle+(X_{2}-v_{})\left(r_{2}(v_{}-X_{2})+a_{21}(X_{1}-u_{*})\right)$
	$\displaystyle+\frac{a_{13}a_{21}}{a_{12}a_{31}}(Y-w_{})\left(a_{13}(X_{1}-u_{})\right)-\frac{\rho a_{13}a_{21}}{a_{12}a_{31}}\frac{(Y-Z)^{2}}{Y^{2}}$
$\displaystyle=$	$\displaystyle-\frac{a_{21}r_{1}}{a_{12}}(X_{1}-u_{})^{2}-r_{2}(X_{2}-v_{})^{2}-\frac{\rho a_{13}a_{21}}{a_{12}a_{31}}\frac{(Y-Z)^{2}}{Y^{2}}\leq 0.$

Global Dynamics and Existence of Traveling Wave Solutions for A Three-Species Models

Abstract

1 Introduction

2 Global Dynamics of the Corresponding ODE System

2.1 Preliminaries

Lemma 2.1.

Proof.

2.2 Assumptions and Two Extinction Results

Lemma 2.2.

Proof.

Lemma 2.3.

Proof.

2.3 Equilibria and Stability in ℝ3\mathbb{R}^{3}

Proposition 2.1.

Theorem 2.1.

Proof.

Theorem 2.2.

Proof.

3 Existence of traveling wave solutions

3.1 The ODE forms and the Lienard Transformation

3.2 The Wazewski set and its Exist Subset

Lemma 3.1.

Proof.

Lemma 3.2.

Proof.

3.3 The Unstable Manifold of E1E_{1}

Lemma 3.3.

Proof.

Lemma 3.4.

Proof.

3.4 Existence and Non-existence of Traveling wave solutions from E1E_{1} to E∗E_{*}

Proposition 3.1.

Proof.

Proposition 3.2.

Proof.

4 Numerical Simulations and Brief Discussions

4.1 Numerical Simulations

4.2 Brief Discussions

References

2.3 Equilibria and Stability in $\mathbb{R}^{3}$

3.3 The Unstable Manifold of $E_{1}$

3.4 Existence and Non-existence of Traveling wave solutions from $E_{1}$ to $E_{*}$