INTRODUCTION
The modeling of preypredator interactions is of great importance in the study
of population dynamics and mathematical ecology, since the early days of ecological
science discipline. The basic preypredator interaction is normally described
by using a system of ordinary differential equations which modeled the spatial
distribution of species as the time evolves (Wang et al.,
2007).
The presence of the diffusion mechanism in preypredator interaction however
changes the behavior and nature of the whole model. It is now a Partial Differential
Equation (PDE) which can be categorized as a reactiondiffusion system. Actually,
the diffusive system has been the focal point of several analytic works. For
instance, Kaipio et al. (1995) have examined
this model for the case of one species and two species in one dimensional diffusion.
They approximated the solutions by using numerical analysis method, namely Galerkin
finite element method and performed simulations on heterogeneity of environments.
The diffusive preypredator model has also been studied extensively by Dunn
et al. (2009), Britton (1986), Murray
(1989), Levin et al. (1993), Malchow
(1993), Medvinsky et al. (2002) and Holmes
et al. (1994).
Apart from that, the inclusion of diffusion terms has also made our preypredator
model tend to be complicated and it becomes a nonlinear system that is very
difficult to analyze and solve analytically. This means that it does not have
closed form solutions and thus many researchers focus their studies on the existence
and uniqueness of solutions in the long run. There are several powerful methods
that can be used to study the existence of solutions of the preypredator model,
namely the method of invariant region and energy estimates. There are numerous
standard literatures on these concepts (Smoller, 1992;
Logan, 2008; Tveito and Winther, 2005).
In general, this study was concerned with another situation. The main objective
was to predict the large time behavior of solutions by using the notions of
invariant region, energy estimates and exponential stability. Plus, this study
also intended to explore and interpret the occurrence of energy decaying exponentially
in the long run from the ecological point of views.
MATHEMATICAL MODEL
In this study, the diffusive preypredator model takes the form:
with initial condition:
and Neumann boundary condition:
The initial and boundary condition represents the evolution of prey and predator
population in the interval [0, 1]. The growth is constrained by the capacity
of the environment which is normalized to 1. The diffusion term u_{xx}
allows the prey population to move from the part of the domain that has high
population density to the ones with lower densities. The homogeneous Neumann
boundary condition is used because we assumed that the domain is closed and
no migration occurs across boundaries.
In a system of reactiondiffusion equations, the notions of invariant region,
energy estimates and exponential stability are in fact an important matter in
predicting the large time behavior of solutions. Thus, this study was concerned
with these subjects in a diffusive preypredator model.
INVARIANT REGION AND ASYMPTOTIC BEHAVIOUR OF SOLUTIONS
The notion of an invariant region (or invariant interval in one space variable)
is an important idea in this study. It provides a suitable theoretical foundation
and framework for studying large time behaviour of solutions (Chueh
et al., 1977).
To begin with, consider the diffusive preypredator model (1) being written
in the vector form:
where, xεI, t>0 and:
Denote
the solution vector of Eq. 1 and
is the vector field of nonlinear reaction terms. The positive constant D is
the diffusion constant and I is an interval in ,
possibly all of . The functions F and G are continuous and based on
Eq. 1:
F (u, v) = 1v, G (u, v) = u1 
The initial conditions are as below:
If I is not all of , then it can be assumed that Dirichlet or Neumann
boundary conditions is imposed on the boundary. Now, consider the following
definition.
Definition 1: Let Σ be a closed set in ^{2}.
If
is a solution to Eq. 4 and 5 for 0≤t≤δ≤∞,
with given boundary conditions, plus the initial and boundary values Σ
are in and
is in Σ for all x where, xεI and 0<t<δ then is called an
invariant region/invariant set for the solution .
Actually, there is a simple condition on Eq. 4 which guarantees
that a particular region is invariant. Such region is invariant if the reaction
points inwards along the boundary of a rectangle. The subsequent theorem asserts
this condition.
Theorem 1: Let Σ = [a, b]x[c, d] be a rectangle in UV space, with
Σ^{0} denote its interior and
denote the boundary, with
the outward unit normal. If:
then Σ is an invariant set for Eq. 45.
Proof: Theorem 1 will be proven by using contradiction. Let suppose
that Σ is not an invariant set. Assume without loss of generality that
u (x_{0}, t_{0}) = b for some (x_{0}, t_{0})
with u(x, t)<b for all xεI and 0≤t<t_{0} and yet:
So, the function (x, t_{0}), regarded as a function of x, must have
a maximum at x = x_{0}. Hence, Δu (x_{0}, t_{0})≤0.
Furthermore, at (x_{0}, t_{0}):
The latter implication follows from the fact that:
on the boundary u = b. But Eq. 8 contradicts Eq.
7 and thus Σ is an invariant set (Q.E.D).
Equivalently, Theorem 1 states that the projection of the reaction on the outward
normal to the surface of the rectangle is nonpositive. Theorem 1 can be applied
to domains containing points where
or where the outward unit normal
is undefined provided that at these points the direction of the derivative vector
does not point out of the invariant set.
Next, let study the following corollary which can be used to obtain global
existence of solutions and thereby provides a suitable framework for studying
the large time behaviour of solutions.
Corollary 1: If the system admits a bounded invariant region Σ
and u_{0}εΣ (u_{0} is the initial data) for all xε ,
then the solution exists for all t>0.
The proofs of Corollary 1, Theorem 1 and Definition 1 follow directly from
Logan (2008).
Now, let show that Eq. 1 admits a bounded invariant region
as t>0. Based on Theorem 1, define:
Draw the zero sets of the components of vector field
(refer to Fig. 1) claim that:
Σ = {(u, v): 0≤u≤1, 0≤v≤1} 
is invariant. Σ is an invariant region if
(by Theorem 1).
Now, begin to prove the abovementioned claim by considering the edge E_{1}
(refer to Fig. 1). E_{1} is the edge u = 0 with 0≤v≤1.
This edge has
or outward pointing unit normal of:
Then:
For edge E_{2 }(E_{2} is the edge 0≤u≤1 with v = 1),
it has
or outward pointing unit normal of:
Then:
The remaining edges E_{3} (E_{3} is the edge u = 1 with 0≤v≤1)
and E_{4} (E_{4} is the edge 0≤u≤1 with v = 0) can be
similarly checked using the above calculation. Observe that
point inward along the boundaries of the rectangle, thus Theorem 1 holds and
Σ is invariant. By Corollary 1, the solution exists for all t>0.
Next, the notion of invariant region will be used to prove the exponential
stability in a diffusive preypredator model (Fig. 1).

Fig. 1: 
Invariant region for preypredator model (1), E_{1}
is the edge u = 0 with v = 1; E_{2} is the edge 0≤u≤1 with
v = 1; E_{3} is the edge u = 1 with 0≤v≤1 and E_{4}
is the edge 0≤u≤1 with v = 0. f, F and G as in Eq. 9 
EXPONENTIAL STABILITY: ENERGY DECAY
Here, the energy decay in a diffusive preypredator model will be shown by
using the concept of energy method, an important technique that is applicable
to all types of PDE. Firstly, the energy integral is established and then this
integral is then manipulated by using integration by parts and a set of key
inequalities (Raposo et al., 2008).
The general preypredator model of population dynamics with added diffusion
of two species is given below:
with initial condition:
and Neumann boundary condition:
D is a diffusion constant and 0≤x≤1, t≥0.
The quantity E(t), defined by:
is called the energy at time t. This can be interpreted as the energy of a
solution as the square of L^{2}norm of solutions u and v.
Now, an important theorem in this study is outlined as below:
Theorem 2: Consider the Eq. 10 and 11
in Ω, with boundary conditions as in Eq. 13. Assume
that Eq. 10 and 11 admits a bounded invariant
region Σ and that
then there exist positive constant C and σ, such that:
Proof: This theorem shall be proven with respect to the Eq.
10 and 11. By multiplying Eq. 10 by
u and integrating over xεΩ, this will result in:
Next, multiply Eq. 11 by v and integrate over xεΩ,
to get:
Notice that the terms:
and:
are obtained by using integration by part. Adding Eq. 16
and 17, will result in:
whereby the Poincare’s inequality is applied so as to get the above result.
Lastly, by Gronwall’s inequality:
This is the end of proof of Theorem 2 (Q.E.D.).
Note that, the decay of energy depends on quantity 2 (1Dπ^{2})
whereby, if:
then E (t) decays exponentially in the long run. In Theorem 2 it is shown the
occurrence of energy decay in the long run. Energy decaying exponentially means
that E→0 as t→∞. Since, the energy is defined as a solution
as the square of L^{2}norm of solutions u and v, this indicates that
both solutions u and v would also decay to zero as t goes to infinity.
The following definition by Korman (1990) explains the
above phenomenon of u and v decay to zero.
Definition 2: It is said that a species u (x, t) dies out or extinct
if
otherwise, u (x, t) persists.
NUMERICAL RESULTS AND THE SIGNIFICANCE OF ENERGY DECAY
The model (1) is discretized using finite difference method whereby u_{j,
m} and v_{j, m} denotes the approximations of u (x_{j},
t_{m}) and v (x_{j}, t_{m}), respectively:
with:
j = 0, 1, …, n + 1 and m≥0. The Neumann boundary condition is as follows:
For initial conditions:
And, choose initial conditions such that:
and v (x, 0) = cos^{2} (5πx).

Fig. 2(ab): 
Extinction of population when D = 0.2, (a) Prey and (b) Predator 
Now, several numerical investigations are conducted in order to corroborate
the theoretical results with numerical evidences. Based on Theorem 2 and Definition
2 the exponential decaying of energy can be interpreted as the extinction of
prey and predator populations in one ecosystem. This occurs when:
Whilst if D is too small i.e.:
the decaying of energy is not guaranteed and a population persistence can occur
in our preypredator model (Raposo et al., 2008).
Consider finitedifference solution plots below with D = 0.2 (extinction) and
D = 0.01 (persistence):
Based on the Fig. 2, observe that the populations of preypredator
have fallen at the end of simulation and settle down to a certain value near
zero. Thus, the extinction of population occur when:
And from Fig. 3, the population of preypredator persist
when:
The phenomena of extinction and persistence of species is also being analyzed
by Gopalsamy (1977) using the model of competing species:

Fig. 3 (ab): 
Persistence of population when D = 0.01, (a) Prey and (b)
Predator 
Gopalsamy analyzed model (20) and obtained the conditions for existence of
nonnegative solutions for all t. For the sake of discussion, let outline the
main result below.
Assume that:
Eq. 20 are subjected to appropriate initial and boundary
conditions at x = 0 and x = L.
It is discovered that the global extinction occur if:
On the other hand, if:
the persistence of species is possible. In particular, the two competing species
can survive and coexist in the ecosystem.
Another interesting remark on the exponential stability is about the large
time behavior of solutions, satisfying Neumann boundary conditions. Theorem
2 indicates that if there is an invariant region Σ and Eq.
15 holds, then all solutions decay as t→∞ to spatially homogeneous
solutions.
In other words, under suitable circumstances, it might be reasonable to ignore
the diffusion process and assume that u and v does not vary too much from point
to point in space. This will lead to system of ordinary differential equations:

Fig. 4(ab): 
Surface plots of model (1) for (a) Prey and Predator population 

Fig. 5(ab): 
Surface plots of model (1) for (a) Prey and (b) Predator
population when D = 0.001 
Then, the solutions to diffusive preypredator model can be approximated by
the abovementioned system of ordinary differential equations. This is sometimes
referred to as the “lumped parameter assumption”. Smoller
(1992) for the discussion on this matter.
In order to see this, consider the solution plots below for the case D = 0
(no diffusion) and D = 0.001 (very small diffusive constants).
Based on Fig. 4 and 5, similar surface
plots of solutions for the case D = 0 (no diffusion) and D = 0.001 (very small
diffusive constants) are obtained. The diffusive preypredator model can be
approximated with system of ordinary differential equations in the long run
when the diffusion mechanism is too small. And this is consistent with the “lumped
parameter assumption”.
CONCLUSION
In this study, the invariant region of the diffusive preypredator model has
been constructed and proven. One of the significant roles of invariant region
is to predict the long run behaviour of solutions of our reactiondiffusion
equation. If such invariant region can be found, then one automatically obtains
a priori bounds on the solution and this priori bounds are essential in obtaining
global existence of solutions. Besides that, the occurrence of energy decaying
exponentially in the long run is proven by using classical energy method. Energy
decaying exponentially means that E→0 as t→∞. Since, the energy
is defined as a solution as the square of L^{2}norm of solutions u
and v, this indicates that both solutions u and v would also decay to zero as
t goes to infinity. This will result in the extinction phenomenon of preypredator
populations; otherwise will result in persistence of population. All these are
very important to be observed and studied in order to understand the ecological
interactions between prey and predator.
ACKNOWLEDGMENT
Research supported by USM Fundamental Research Grant Scheme (FRGS).