Spatial Patterns of a Predator-Prey System of Leslie Type with Time Delay.

Time delay due to maturation time, capturing time or other reasons widely exists in biological systems. In this paper, a predator-prey system of Leslie type with diffusion and time delay is studied based on mathematical analysis and numerical simulations. Conditions for both delay induced and diffusion induced Turing instability are obtained by using bifurcation theory. Furthermore, a series of numerical simulations are performed to illustrate the spatial patterns, which reveal the information of density changes of both prey and predator populations. The obtained results show that the interaction between diffusion and time delay may give rise to rich dynamics in ecosystems.

Introduction

Thanks to the classical work of Lotka (in 1925) and Volterra (in 1926), modeling predator-prey interaction system has become one of the hot issues in mathematical ecology [1-6]. As is well known, one of the principles that predator-prey models follow is that predators can grow as a function of what they have eaten [7]. One of the famous functional response function is Generalized Holling type III [8]. When b = 0, it is called Holling type III. Moreover, not only is the predator growth term described by a function of the prey density, but also is described as a function of the ratio of predator and their prey, y/x, where x and y stand for prey and predator density respectively, see for example [9, 10]. The predator-prey system takes the following form: where g(x, K) describes the specific rate of the prey if there is no predator. p(x) is the functional response function which describes the change in the density of the prey when they are attacked by per predator in per unit time. Information about the properties of function g(x, K), p(x) and q(x) are available in [11, 12]. In this paper, we consider system (1) with the following functions, and , namely, the predator-prey system of Lesile type with Holling type III functional response. Then Eq (1) becomes: where the parameter r, K, a, θ, m and h are all positive constants. r is the prey intrinsic growth rate. K is the carrying capacity. m is capturing rate. a is half capturing rate. θ is predator intrinsic growth rate. h is conversion rate of prey into predator biomass. For the sake of convenience, system (2) should be rewritten into the nondimensional form. Assuming then system (2) becomes:

Pattern formation in reaction-diffusion system is one of the attractive problems in natural, social, and technological sciences. Espatially in ecological system, various predator-prey models with diffusion have been studied [13-18]. Pattern formation can well explain species survival under the influence of individual mobility. Combing Eq (3) with diffusion, we have the spatiotemporal predator-prey system of Lesile type with Holling type III functional response: where the positive constants D1 and D2 denote the diffusive coefficients of u and v, respectively. is the usual Laplacian operator in two-dimensional space, which describes the random motion.

On the other hand, time delay due to maturation time, capturing time, gestation or other reasons widely exists and plays an important role in many biological dynamical systems [19, 20]. In order to reflect the current population dynamics, the rate of change of which depends on the past population of the system, we should incorporate time delays into mathematical models [21-24]. Although a lot of work has been done about the spatial predator-prey model [25-28] and studies of delay feedback on pattern formation have achieved great progress [29-33], study of delay driven pattern formation in a Leslie type system with Holling type III functional response seems to be rare. As a result, in the present paper we aim to study the effects of time delay on the spatiotemporal dynamics of a Leslie type model with Holling type III functional response. The model is as following form: where τ > 0 is a constant due to the negative feedback. Let Ω be a square flat domain. The initial conditions are Generally speaking, to make sure that Turing pattern is determined by reaction-diffusion mechanism, we usually choose zero-flux boundary conditions which means that there is no flux of populations through the boundary, i.e., no external input is imposed from outside.

This paper is organized as follows. In section 2, we study the dynamics of model without delay. In section 3, we obtain the condition of Turing instability (diffusion and delay induced instabilities)via linear stability analysis. In section 4, we present various spatial patterns by performing numerical simulations. Finally, we give some conclusions and discussions in section 5.

Materials and Methods

Existence of positive equilibria

We need to analyze the stability criteria of model (5) without delay and diffusion. The corresponding model is system (3). Obviously, system (3) has equilibrium E0 = (1, 0), which corresponds to extinction of the predator. From the biological point of view, we are interested in the interior equilibria points, which are the positive solutions of the following cubic polynomial equations of the system (3):

Substituting the second equality for the first equality in Eq (8), we have that: where and . The number of equilibria in Eq (3) is determined by the number of real roots of F in the interval I0 = (0, 1). In addition, F′(u) = 3u2 + 2ω2u + ω1 has two zeros when The discriminant of the cubic polynomial F is given by where and . Similar to Lemma 1 in [34], we have the following results.

System (3) has at least one equilibrium and at most three equilibria in the interval I0 = (0, 1). Moreover, consider the following condition:

If and only if Eqs (13) and (14) hold, then system (3) has three equilibria in the interval I0 = (0, 1).

If and only if Eqs (15) and (16) hold, system (3) two equilibria in the interval I0 = (0, 1).

If and only if Eqs (18), or(13) and (17) hold, system (3) has a unique equilibrium in the interval I0 = (0, 1) Moreover, if system (3) has a unique equilibrium where , then is either a simple zero of F or a zero of multiplicity 3 of F.

We describe qualitative properties and stability of the interior equilibria of system (3) (see S1 File). By the proof of Theorem 1 (see S1 File), we have the following conclusion.

If all of the positive equilibria of system (3) are stable, except the saddle (corresponds to the zero A of F in Fig 1) and degenerate equilibria (corresponds to the zeros B, C and D of F in Figs 2 and 3 respectively).

Fig 1

F has three zero solutions in I0.

Fig 2

F has two zero solutions in I0.

Fig 3

F has a unique zero solution in I0.

Turing instability analysis of model (4)

In this subsection, we will give a brief analysis of Turing instability of model (4). For the sake of convenience, let E*(u*, v*) be anyone of the interior steady equilibrium in Theorem 1 (see S1 File). The characteristic polynomial at E*(u*, v*) is where J = J − diag(D1, D2)k2, K is a wavenumber and J is the Jacobian matrix of system (3) at E*(u*, v*).

Further, Eq (20) yields where and

We can get the roots of Eq (20):

The condition for the onset of Hopf instability holds when a pair of imaginary eigenvalues cross the real axis from the negative value to the positive one and there is no diffusion [35, 36]. That is to say, the Hopf bifurcation occurs when

A homogeneously steady state is said to be Turing instability if it is stable for model (3) without diffusion but becomes unstable because of homogeneous perturbation caused by diffusion. A general linear analysis [37-39] show that the necessary conditions for onset of Turing instability for model (4) are given by condition Eq (19) and the following conditions The condition Eqs (19) and (25) make sure that the equilibrium E* = (u*, v*) is stable for model (3) without diffusion, and becomes unstable for model (4) if Re(λ) transits the real axis from a negative side to a positive one (corresponding to condition Eq (26) and (27)). Namely, the Turing bifurcation occurs when and wavenumber k satisfies

Results

Linear stability analysis

In this subsection, we will consider the stability of model (5). Obviously, model (5) has the same equilibria as model (4). Similar to [30, 31], assume that τ is small enough, then we replace and in system (5) and obtain the following equations: Expanding Eq (28) in Taylor Series and neglecting the higher order non-linearities, then Eq (28) becomes: where and . By Eq (29), we obtain the following equations:

If we take small spatiotemporal perturbations δu(x, y, t) and δv(x, y, t) on the steady state E* = (u*, v*) of system (5), then we have:

Expanding the reaction terms around the steady state E* = (u*, v*) in Taylor Series up to first order and rearranging the terms, we obtain: where with . Since τ is small, we only consider in this paper. Model (32) becomes:

Assume that spatiotemporal perturbations δu(x, y, t) and δv(x, y, t) take the following form: where λ is the growth rate of the perturbation in time t, δu* and δv* stand for the amplitudes, and k and k are the wavenumbers of the solutions. Inserting Eqs (34) into (33), we obtain the characteristic equation at E* = (u*, v*) of model (5): where

Now we are interested in investigating the effects of time delay and diffusion on the dynamical system (5), and we want to know under what conditions for time delay to destabilize the steady state and let the spatiotemporal instability occur. The onset of Turing instability requires at least one of tr(J) < 0 and det(J) > 0 is violated. So we consider the emergence of the Turing instability in the following two cases:

det(J) > 0 is violated.

tr(J) < 0 is violated.

Diffusion induced instability

In this part, we will consider the first case, i.e. det(J) > 0 is violated, namely, det(J) < 0. It can be seen from Eq (38) that det(J) = χ det(J) and , so the sign of det(J) is the same as det(J). From Eq (23), det(J) < 0 is equivalent to: at the critical value of wavenumber .

In addition, tr(J) < 0 equals to . Because that , simple algebraic computation leads to: Hence, in this case, τ must satisfy

Finally, we can get the condition of diffusion induced instability:

When system (3) has three equilibria, if condition Eqs (11), (13), (14), (19), (39), (40), and (42) hold, then instability of model (5) induced by diffusion occurs.

When system (3) has two equilibria, if condition Eqs (11), (15), (16), (19), (39), (40) and (42) hold, then instability of model (5) induced by diffusion occurs.

When system (3) has a unique equilibrium, if condition Eqs (11), (18) (or (13) and (17)), (19), (39), (40) and (42) hold, then instability of model (5) induced by diffusion occurs.

Delay induced instability

In this part, we will consider the seconde case, i.e. tr(J) < 0 is violated. In order to find suitable values of parameters for our simulation, we assume in this part. Following the same analysis in subsection B, tr(J) < 0 is violated when , which yields Hence, in this case, τ must satisfy:

We mainly consider the delay induced instability, therefore we keep det(J) > 0, which is equivalent to the same condition as Eqs (39) and (40).

At last, we can get the condition of delay induced instability:

When system (3) has three equilibria, if condition Eqs (11), (13), (14), (19), (39), (40) and (44) hold, then instability of model (5) induced by delay occurs.

When system (3) has two equilibria, if condition Eqs (11), (15), (16), (19), (39), (40) and (44) hold, then instability of model (5) induced by delay occurs.

When system (3) has a unique equilibrium, if condition Eqs (11), (18) (or (13) and (17)), (19), (39), (40) and (44) hold, then instability of model (5) induced by delay occurs.

In this section, we will perform numerical simulations for model (5) on the 100 × 100 square lattices with Neumann boundary conditions. The simulations are initiated with small amplitude random perturbations around the positive equilibrium point E(u*, v*). The reaction-diffusion equations in our models are analyzed numerically employing Forward Difference implicit difference scheme. We set the time step Δt = 0.005, spatial mesh size h = 1.

We run the simulations until they reach a stationary state which indicates that the behavior does not seem to change its characteristic anymore. As a result, we only make analysis of pattern formation to one distribution(in this paper, we show the distribution of the predator).

Numerical results

A. Pattern formation induced by diffusion

In this subsection, we will focus on the pattern formation induced by diffusion, we obtain nine different types of patterns in Fig 4, where prey and predator can coexist.

Fig 4

Nine types of patterns induced by diffusion.

In Fig 4(a) we take the parameters as Fig 4(a) is the typical spot pattern, which shows that the distribution of the predator v is the isolated regions with high density.

In Fig 4(b) we take the parameters as Fig 4(b) is the typical hole pattern, which shows that the distribution of the predator v is the isolated regions with low density.

In Fig 4(c) we take the parameters as Fig 4(c) is the mixture of red stripes and red spots named as mixed pattern 1.

In Fig 4(d) we take the parameters as Fig 4(d) is the mixture of blue stripes and blue spots named as mixed pattern 2.

In Fig 4(e) we take the parameters as Fig 4(e) is the stripe pattern.

In Fig 4(f) we take the parameters as Fig 4(f) is somewhat like labyrinth pattern.

In Fig 4(g) we take the parameters as Fig 4(g) is somewhat like grid pattern.

In Fig 4(h) we take the parameters as Fig 4(h) is the composed parallel lines.

In Fig 4(i) we take the parameters as Fig 4(i) is somewhat butterfly-like patten.

B. Pattern formation induced by delay

In this subsection, we will focus on the pattern formation induced by delay, we obtain three different types of patterns in Fig 5.

Fig 5

Three types of spirals induced by delay.

In Fig 5(a) we take the parameters as Fig 5(a) is spirals 1.

In Fig 5(b) we take the parameters as Fig 5(b) is spirals 2.

In Fig 5(c) we take the parameters as Fig 5(c) is spirals 3.

Discussion

In this paper, a spatial Lesile type predator-prey system with Holling type III functional response and time delay has been investigated. To well understand the impact of delay and diffusion on the instability, we have made theoretical analysis and numerical simulations. Since the equilibrium cannot be expressed in a useful closed form, we cannot discuss its qualitative properties in normal routine. Firstly, we discuss the number and qualitative properties of positive equilibrium via the original parameters. Secondly, we obtain conditions of two types of instability: diffusion induced instability and delay induced instability. Finally, numerical simulations are performed to illustrate the theoretical findings. Both the theoretical and numerical results reveal that the interaction between time delay and diffusion can give rise to stationary patterns.

However, it should be noted that the method in this paper is only suitable for short time delay τ. When the delay is large, one should use other methods to find the condition for Turing instability. Moreover, we may investigate travelling wave of model (5) in the future study. It should be also worth pointing that other types of instability may be found in system (5). For example, we can use normal formal theory and the center manifold theorem of partial functional differential equations to analyze the Hopf bifurcation of system (5) [3].

Qualitative properties and stability of the interior equilibria of system (3).

(PDF)

Click here for additional data file.