Analysis and simulation on dynamics of a partial differential system with nonlinear functional responses

Abstract. We introduce a reaction–diffusion system with modified nonlinear functional responses. We first discuss the large-time behavior of positive solutions for the system. And then, for the corresponding steady-state system, we are concerned with the priori estimate, the existence of the nonconstant positive solutions as well as the bifurcations emitting from the positive constant equilibrium solution. Finally, we present some numerical examples to test the theoretical and computational analysis results. Meanwhile, we depict the trajectory graphs and spatiotemporal patterns to simulate the dynamics for the system. The numerical computations and simulated graphs imply that the available food resource for consumer is very likely not single.


Introduction
Many ecological phenomena among different populations can be characterized or simulated by various mathematical models. By analyzing different kinds of mathematical models people may give scientific predictions and explanations on the dynamics of these models, and further, put forward reasonable projects corresponding to some ecological problems. In many kinds of ecological-mathematical models, the predator-prey model is a very important branch. Since the classical Lotka-Volterra ecological model [13] was introduced into investigations, the predator-prey-type models received extensively attentions. In recent decades, different kinds of predator-prey models have been established, and many great progress and valuable achievements have been made; see [4,6,9,20,22,23] for examples.
In predator-prey models, the functional responses are very important terms in describing the relations between predator and prey, they determine many dynamical properties of the systems in some extent. Specifically, the functional response describes the change of predation rate of predator capturing prey, that is, how the predator consuming prey is dependent on the response terms. In large number of predator-prey models, the functional responses are of prey-dependent type, such as the classical Lotka-Volterra-type functional response, the Holling-type II (sometimes referred to as the Michaelis-Menten functional response), III and IV functional responses, the Ivlev and inverse Ivlev functional responses, etc. Taking these functions as the functional responses, a great many of predator-prey models have been studied thoroughly, and quite a lot of important research achievements have been given; see [1,10,18,21,24], etc.
The introduction of prey-dependent-type functional response is based on the assumption that the predation rate of predator is only effected by prey, and the predators are mutual noninterference within themselves. In fact, generally speaking, the phenomenon in sharing or competing their foods always occurs during the process of predation. When the population density of predator is high, the predation rate will decrease. Moreover, as the well-known environmental paradox being proposed by Hairston, et al. [8] and Rosenzweig [19], the type of functional responses gave rise to many strong controversies among ecologists. Just in this situation, based on the Holling-type II functional response, by a great deal of experimental and numerical analysis, Beddington [3] and DeAngelis et al. [5] introduced a new type functional response, the predator-prey dependent functional response f (x) = bx/(a + x + cy), which is called the Beddington-DeAngelis-type functional response, where x and y represent the predator and prey, respectively, and a, b, c are parameters. This type of functional response is very like with the Hollingtype II functional response, the only difference is that the Beddington-DeAngelis-type functional response involves both prey and predator. Compared with the prey-dependenttype functional responses, the predator-prey-dependent functional response not only reflects the mutual interference with predators, but also remains the merits of the ratiodependent functional responses. Moreover, it avoids some controversies induced by low population density and better describes the predation effect of the predator versus prey. It is just because of these reasons, different kinds of models with predator-prey-dependent functional response have become a center topic in population dynamics.
Except the Beddington-DeAngelis-type functional response, another attentive predator-prey-dependent functional response hx/y is called the Leslie-Gower-type functional response, where h is a constant. It was introduced by Leslie [11] and discussed jointly by Leslie and Gower [12]. This functional response measures the self-consumption of predator when the foods are scarce, where the growth of predator population is of logistic form. Since in actual ecological systems, the carrying capacity set by environmental resources is proportional to prey abundance and the predators always try to survive by catching other preys when their conventional food is short seriously, in such situations, the growth rate of the predator would be affected. For this reason, Aziz-Alaoui [2] modified the logistic form and proposed the Leslie-Gower-type functional response hx/(a + y), where a is a positive constant and measures the environmental protection for predator. With such a functional response, Aziz-Alaoui et al. studied some predator-prey models in spatially homogeneous or inhomogeneous cases; see [7,[15][16][17]25]. Based on and motivated by all the above-mentioned, in this paper, we introduce the following nondimensional reaction-diffusion system with two modified predator-prey-dependent nonlinear functional responses αv/(c + u + mv), βv/(r + u).
Here Ω is a bounded domain in R n with smooth boundary ∂Ω, u and v are the densities of prey and predator populations. d 1 , d 2 and a, b are the diffusion rates and the growth rates of u and v, respectively; α is the predation rate of predator, and c is the self-saturation density of prey; m represents the semi-saturation term of the system; β measures the conversion rate of v predating u, and r reflects the protection of the surroundings on predators; ν is the outward unit normal vector on ∂Ω. All parameters are positive constants. u 0 (x), v 0 (x) are booth smooth functions in Ω.
For system (1), we consider the large-time behavior of positive solutions. In addition, corresponding to (1), we mainly analyze the steady-state system An outline of this paper is as follows: In Section 2, we consider the large-time behavior of positive solutions to system (1) mainly by the comparison principle of parabolic equations. In Section 3, we give a priori estimate for the positive solutions of (2) by the maximum principle. Then, in Section 4, we investigate the coexistence of predator and prey to system (2) by the energy integral techniques and topological degree computation. In Section 5, by taking the diffusion rate of predator as a parameter, we study the bifurcating phenomenon, which emits from the unique positive constant equilibrium also by using the topological degree techniques. Since pattern is a very interesting nonlinear phenomena and pattern dynamics is a primary branch of nonlinear science, finally, we give some numerical examples and depict the corresponding trajectory graphs or spatiotemporal patterns to simulate the related theoretical results in Section 6. It is worth mentioning that the numerical simulated graphs imply that the available food resource for predator may not single if the system parameters are controlled properly.

Large-time behavior
In this section, we consider the large-time behaviors of positive solutions to system (1) including the global attractor and the persistence.
with Ω being the same as that of in (1) and f 0 being a smooth function in Ω, then lim t→∞ f (x, t) = k holds uniformly in Ω.
Likewise, using the comparison principle of parabolic equations and Lemma 1 again, we also know that there is Then the arbitrariness of ε induces that the conclusion holds.
Proof. For any ε > 0, the proof of Theorem 1 implies that u(x, t) satisfies where t 2 is defined in the proof of Theorem 1. Similar as that we prove Theorem 1, there Then the arbitrariness of ε implies that This shows that system (1) is persistent.
Remark 1. Theorem 2 implies that the predator and prey always coexist in their living surroundings despite the location and time as long as ac > bα, moreover, this coexistence is independent of their diffusion situation.

A priori estimate
In this section, we seek a priori estimate for the positive solutions of (2), which will be used frequently in the latter contents. As preparation, we introduce the following lemma due to [14].
Then the followings hold.
, v(x)) be any positive solution of (2). Then Proof. Since (u(x), v(x)) satisfies (2), then regularity theory of elliptic equations implies that u(x) and v(x) must attain their maximum and minimum values in Ω. By the first equation of (2) we have Then Lemma 2(i) induces u(a − u) 0 at the maximum of u, so maxΩ u(x) a. Then by maxΩ u(x) a and the second equation of (2) we get Using Lemma 2(i) again, there holds v(b − v − βv/(a + r)) 0 at the maximum of v, and further, maxΩ v(x) (a + r)/(a + r + β) < b. The proof is finished.
Similar to the proof of Theorem 3, apply Lemma 2(ii) directly to the first equation of (2) to yield Thus, where P > 0 is a constant satisfying P > bα(r + u(x 2 ))/(cu(x 2 )(r +β +u(x 2 ))). Now, let Then h(x) and u satisfy Further, by the Harnack inequality when Substitute (6) into (5) to get Combining (4) with (7), we have The result follows.
In a same way, if we take d 2 as the parameter, then we can get a very similar estimate on the lower bound.

Nonexistence
This subsection devotes to the nonexistence of nonconstant positive solutions to system (2) by energy integral procedure. The following two nonexistent results are based on the priori estimates obtained in the previous section. Proof. We only prove the result holds under the first condition since the second case can be proved similarly. Suppose that (u, v) is a nonconstant positive solution to system (2). Integrating the first equation of (2) in Ω, we get Using the positivity of u, we obtain further in Ω to get Similarly, we can get Add these two inequalities and use the Cauchy inequality to yield where ε > 0 is arbitrary and L = max{α, βv 2 /((r + C)(r + u))}. Using the Poincaré equality, we have Since d 2 > bµ −1 1 , take ε small enough such that µ 1 d 2 > b + εL. Then there holds Then this inequality holds if and only if u ≡ u is a constant in view of d 1 > D. This is a contradiction, and the proof is complete.

Existence
This subsection deals with the existence of nonconstant positive solutions to system (2) by choosing d 2 as the parameter and by using the Leray-Schauder degree theory.
Firstly, we give some preliminaries. Denote Denote by 0 = µ 0 < µ 1 < µ 2 < · · · the eigenvalues of −∆ under homogenous Neumann boundary condition in Ω and E(µ i ) the eigenspace corresponding to each with So, U is a solution of (2) if and only if U is a solution of (8), this is equivalent to In the following, we need to use the positive constant solution of (2). By a direct analysis on the derivative of a cubic function it can be easily shown that (2) admits a unique positive constant solution denoted by (u * , v * ) =: U * if one of the following two cases holds.
Clearly, U * is a zero point of G(d 1 , d 2 ; ·), and the Fréchet derivative of G(d 1 , It easy to check that λ is an eigenvalue of D U G(d 1 , d 2 ; U * ) on X i if and only if λ(1 + µ i ) is an eigenvalue of the matrix X il it is seen that λ is an eigenvalue of D U G(d 1 , d 2 ; U * ) on each X il , also, the multiplicity of λ is equivalent to the multiplicity of the eigenvalues λ(1 + µ i ) of M i . Therefore, if λ is an q-multiplicity eigenvalue of D U G(d 1 , d 2 ; U * ) on each X i and λ(1 + µ i ) is an n-multiplicity eigenvalue of M i , then q = nm(µ i ).
These contradict to (17), and the proof is completed.
and then it is equivalent to W satisfying The Fréchet derivative of R(d 2 ; ·) on W at W = 0 is where B is given in Section 4.2. As what we do in Section 4, λ is an eigenvalue of and M i is nonsingular, so D W R(d 2 ; 0) is regular. By the implicit function theorem we know that there exists ε, 0 < ε 1, such that W = 0 is the unique solution of R(d 2 ; ·) = 0 in Ξ = {U ∈ X: U X < ε} for any d 2 (d 2 is in a neighborhoodd 2 ), this means that (d 2 ; U * ) is a regular point of (2).
(ii) By assumptions we take µ i ∈ S ∩ N (d 2 ), then H(d 1 ,d 2 ; µ i ) = 0 and µ i = (a 11 d −1 Thus, the two eigenvalues of M i are 0 and tr M i , respectively, and 0 is simple. If the result is false, then there isd 2 > 0 satisfying the following.

Numerical examples, discussions and conclusions
In the previous sections, we investigate the coexistence of prey and predator and find out their coexistent conditions. Nevertheless, from the realistic point of view, we do hope the prey and predator can coexist since our aim is to maintain the balance of ecosystems. The conditions in Theorems 7 and 8 ensure that the prey and predator can coexist and, in this case, the ecological balance may be sustainable. To clarify this point, in this section, we give some numerical examples to support our theoretical analysis. We fix a, b, c, m, r and take α, β as parameters. For the sake of simple calculation, the parameters a, b, c, m and r are referred to [18] in part and taken as follows: a = 1.2, b = 0.08, c = 2.4, m = 50, r = 12. The parameters α, β are chosen, and the values (u * , v * ), ac, bα, brα/ (cr + cβ + bmr) =: s 1 , c + r + β + bm =: s 2 , a 11 are calculated and listed in Table 1. Summarizing these numerical simulations, the following results follow: In Examples 1-2, we see that ac > bα holds, then Theorem 2 shows that system (1) has persistent property and any positive solution (u(x, t), v(x, t)) of (1) satisfies lim inf respectively. However, ac < bα holds in Examples 3-5, so, we do not know whether system (1) has persistence or not in these three cases. Nevertheless, we also see that the existent condition (H1) of (u * , v * ) in Examples 1-5 holds and (u * , v * ) is worked out in each example, and then (2)            All the trajectory graphs show that the quantity of prey and predator decrease and increase gradually as the predation rate of predator, the conversion rate of predator due to capturing prey increasing and time going on, and tend to stabilize eventually. Specifically, the quantity of prey and predator are as followings in turn: 0.08, 0.07, 0.06, 0.05, 0.04 and 0.22, 0.41, 0.65, 0.84, 0.98, respectively, from Figs. 1-5. We also see that the declining and ascending scales between the prey and predator are inconformity, in contrast, the increase of the predator is much stronger than the decrease of the prey. Though the quantity of the prey is decreasing, it is clear that such changes are very small, that is, the quantity of the prey remains roughly stable regardless some fluctuations of the predation rate and the conversion rate of predator. This situation probably implies that the available food resource for predator is not single, that is, what we want to see in population dynamics since which sustains the persistence of the species when the food resource for predator is scarce. This helps to maintain the ecological balance. Figures 6-10 show that the spatiotemporal patterns always occur for our given parameter values, which imply that the system exhibits rich dynamical behavior.   In fact, we made a lot of numerical examples and trajectory graphs, we found that the differences among these simulation and graphic results are very small. Here, we only  present the above five numerical examples and their corresponding trajectory graphs, the spatiotemporal pattern formation to state our findings.
The numerical examples tell us that the ecosystem reflected by model (1) or (2) is easy to maintain stability when the predation rate and the conversion rate of predator do not change dramatically. This might imply that neither the prey nor the predator will disappear within a certain period of time, which could also mean that other species is difficult to invade such system, this is exactly what we hope to happen biologically.