Pattern Formation in a Semi-Ratio-Dependent Predator-Prey System with Diffusion

Abstract

We investigate spatiotemporal dynamics of a semi-ratio-dependent predator-prey system with reaction-diffusion and zero-flux boundary. We obtain the conditions for Hopf, Turing, and wave bifurcations of the system in a spatial domain by making use of the linear stability analysis and the bifurcation analysis. In addition, for an initial condition which is a small amplitude random perturbation around the steady state, we classify spatial pattern formations of the system by using numerical simulations. The results of numerical simulations unveil that there are various spatiotemporal patterns including typical Turing patterns such as spotted, spot-stripelike mixtures and stripelike patterns thanks to the Turing instability, that an oscillatory wave pattern can be emerged due to the Hopf and wave instability, and that cooperations of Turing and Hopf instabilities can cause occurrence of spiral patterns instead of typical Turing patterns. Finally, we discuss spatiotemporal dynamics of the system for several different asymmetric initial conditions via numerical simulations.

1. Introduction

In recent years, pattern formations in nonlinear complex systems have been one of the central problems of the natural, social, technological sciences and ecological systems [1–7]. Particularly, many researchers have studied the prey-predator system with reaction-diffusion. For example, in [8], Garvie and Trenchea presented the analysis of reaction-diffusion predator-prey systems with the Holling type II functional response and provided an 𝐿^∞ a priori estimate. Wang et al. in [9] investigated the spatial pattern formation of a predator-prey system with prey-dependent functional response of Ivlev type and reaction-diffusion. Also Zhang et al. in [10] studied a linear stability and bifurcation analysis including the Hopf and Turing bifurcation for a spatial Holling-type IV predator-prey model. In addition, Camara and Aziz-Alaoui, the authors of [11], considered a predator-prey system with a modified Leslie-Gower functional response modeled by a reaction-diffusion equation and derived the conditions for Hopf and Turing bifurcation in the spatial domain.

In this context, in this paper, we will focus on the following a semi-ratio-dependent predator-prey system with reaction-diffusion: 𝜕𝑁𝑁𝜕𝑇=𝑟𝑁1−𝐾−𝑎𝑁𝑁+𝑏𝑃𝑃+𝐷₁∇²𝑁inΩ,𝜕𝑃𝜕𝑇=𝑒𝑃1−ℎ𝑃𝑁+𝐷₂∇²𝑃inΩ, (1) where 𝑁,  𝑃 stand for the population densities of prey and predator with time 𝑇, respectively. The parameters 𝑟,  𝐾,  𝑎, and 𝑏 are all positive constants standing for the prey intrinsic growth rate, the carrying capacity of prey, the capturing rate, and the half-saturation constant, respectively. The predator grows logically with the growth rate 𝑒 and the carrying capacity 𝑁/ℎ. In fact, the parameter ℎ is a measure of the food quality that the prey provides for conversion into the predator. And 𝐷₁ and 𝐷₂ are diffusion coefficients. Here ∇²=𝜕²/𝜕𝑥²+𝜕²/𝜕𝑦² is the usual Laplacian in two-dimensional and 𝑥 and 𝑦 stand for the space.

In fact, the authors in [12–14] have explored a similar model to system (1) without reaction-diffusion, which is a kind of a temporal predator-prey system, and figured out dynamic properties of the temporal system such as stabilities of equilibria and permanence of the system. However, as mentioned above, it is important to investigate spatiotemporal dynamical behaviors of a diffusive predator-prey system like system (1). Thus, in this paper, we will concentrate our concerns on studying the bifurcation analysis and the spatiotemporal pattern formation analysis of system (1).

Throughout this paper, we assume that no external input is imposed from the outside. Hence the boundary conditions are taken as 𝜕𝑁/𝜕𝑛=𝜕𝑃/𝜕𝑛=0 on 𝜕Ω, Ω is the spatial domain in ℝ², and 𝑛 is the outward unit normal vector of the boundary 𝜕Ω. The main reason for choosing such boundary conditions is that we are interested in the self-organization of pattern. Also system (3) needs to be analyzed with the initial populations 𝑁(𝑥,𝑦,0)>0 and 𝑃(𝑥,𝑦,0)>0.

In order to minimize the number of parameters in system (1), we set 𝑁=𝐾𝑁,𝑃=𝑟𝐾𝑎𝑃,𝑇=𝑇𝑟. (2) For simplicity we will omit the bar notation in 𝑁, 𝑃, 𝑇 in the rest of the paper. Thus we can have the following system containing dimensionless quantities: 𝜕𝑁()−𝜕𝑇=𝑁1−𝑁𝑁𝑃𝑁+𝛼𝑃+𝑑₁∇²𝑁,𝜕𝑃𝑃𝜕𝑇=𝛿𝑃1−𝛽𝑁+𝑑₂∇²𝑃, (3) where 𝑟𝑏/𝑎=𝛼,  𝑒/𝑟=𝛿,  𝑟ℎ/𝑎=𝛽,  𝐷₁/𝑟=𝑑₁,  and 𝐷₂/𝑟=𝑑₂.

For convenience, we set 𝑁(1−𝑁)−𝑁𝑃/(𝑁+𝛼𝑃)≡𝐹(𝑁,𝑃) and 𝛿𝑃(1−𝛽(𝑃/𝑁))≡𝐺(𝑁,𝑃).

The main object of this paper is to look into the spatial dynamic behaviors of system (3). For this, in Section 2, we investigate bifurcation phenomena, Hopf, Turing, and wave bifurcations. In particular, we find out sufficient conditions for the Hopf bifurcation, the Turing instability, and the wave bifurcation of system (3). In Section 3, we give numerical simulations of system (3) with random perturbation initial conditions of the stationary solution of the spatially homogeneous system to investigate various spatiotemporal pattern formations such as spotted, stripelike, and spiral patterns. Finally, we discuss dynamical behaviors of system (3) with several different initial conditions in Section 4.

2. Bifurcation Analysis

In order to investigate pattern formations of system (3), first we must think over the nonspatial system of system (3). In fact, the nonspatial system for system (3) has three equilibria, which correspond to spatially homogeneous equilibria of system (3), in the positive quadrant as follows:(i)(0,0) (total extinct): a saddle point,(ii)(1,0) (extinct of the predator): a saddle point,(iii)a nontrivial positive stationary state (𝑁^∗,𝑃^∗) (coexistence of prey and predator) if 𝛼+𝛽>1, where 𝑁^∗1=1−𝛼+𝛽,𝑃^∗=1𝛽𝑁^∗=1𝛽11−.𝛼+𝛽 (4)From the biological point of view, in this paper, we mainly focus on the dynamics of nontrivial stationary state (𝑁^∗,𝑃^∗).

Now, to perform a linear stability analysis for the nontrivial stationary state (𝑁^∗,𝑃^∗), we must linearize the dynamical system (3) around the spatially homogeneous fixed point (𝑁^∗,𝑃^∗) for small space- and time-dependent fluctuations and expand them in Fourier space. For this, let 𝑁^→𝑥,𝑡∼𝑁^∗𝑒^𝜆𝑡𝑒^{𝑖→𝑘⋅→𝑥},𝑃^→𝑥,𝑡∼𝑃^∗𝑒^𝜆𝑡𝑒^{𝑖→𝑘⋅→𝑥}, (5) where ^→𝑥=(𝑥,𝑦) and ^→𝑘⋅^→𝑘=𝑘²; 𝑘 and 𝜆 are the wavenumber vector and frequency, respectively. Then we can obtain the corresponding characteristic equation as follows: 𝐴_𝑘−𝜆𝐼=0, (6) where 𝐴_𝑘=𝐴−𝑘²𝐷, 𝐼 is the identity matrix and 𝐷=diag(𝑑₁,𝑑₂) is the diffusion matrix and 𝐴 is given by 𝐴=𝜕𝐹𝜕𝑁𝜕𝐹𝜕𝑃𝜕𝐺𝜕𝑁𝜕𝐺𝜕𝑃_(𝑁^∗_,𝑃^∗₎=1−2𝑁^∗−𝛼𝑃^∗2(𝑁^∗+𝛼𝑃^∗)²−𝑁^∗2(𝑁^∗+𝛼𝑃^∗)²𝑃𝛿𝛽^∗2𝑁^∗2𝑃𝛿−2𝛿𝛽^∗𝑁^∗≡𝐹_𝑁𝐹_𝑃𝐺_𝑁𝐺_𝑃. (7) Now, (6) can be solved, yielding the so-called characteristic polynomial: 𝜆²−tr𝐴_𝑘𝜆+det𝐴_𝑘=0, (8) where tr𝐴_𝑘=𝐹_𝑁+𝐺_𝑃−𝑘²𝑑₁+𝑑₂=−𝑑₁+𝑑₂𝑘²1−1−𝛿++𝛽𝛼+𝛽𝛼+𝛽²,det𝐴_𝑘=𝐹_𝑁𝐺_𝑃−𝐹_𝑃𝐺_𝑁−𝑘²𝐹_𝑁𝑑₂+𝐺_𝑃𝑑₁+𝑘⁴𝑑₁𝑑₂=𝑑₁𝑑₂𝑘⁴+𝑑₂+𝑑₁𝑑𝛿−²−𝑑𝛼+𝛽²𝛽𝛼+𝛽²𝑘²𝛿+𝛿−.𝛼+𝛽 (9) The solutions of (8) yield the dispersion relation 𝜆^±_𝑘=12tr𝐴_𝑘±tr𝐴_𝑘²−4det𝐴_𝑘. (10) The reaction-diffusion systems have led to the characterization of three basic types of symmetry-breaking bifurcations—Hopf, Turing, and wave bifurcation, which are responsible for the emergence of spatiotemporal patterns [3, 8, 9, 15–30].

2.1. Hopf Bifurcation

Hopf bifurcation is an instability induced by the transformation of the stability of a focus. In fact, the space-independent Hopf bifurcation breaks temporal symmetry of a system and gives rise to oscillations that are uniform in space and periodic in time. Mathematically speaking, the Hopf bifurcation occurs when Im𝜆_𝑘≠0,Re𝜆_𝑘=0at𝑘=0. (11) Thus we can figure out that the critical value of Hopf bifurcation parameter 𝛿 is equal to 𝛿_𝐻=𝛽𝛼+𝛽²+1𝛼+𝛽−1. (12) We assume that the value 𝛿_𝐻 is greater than 0 for the existence of Hopf bifurcation. In other words, from now on, we will consider the parameters 𝛼 and 𝛽 which are satisfied with the following condition: 𝛼²+2𝛽−1𝛼+𝛽𝛽−2<0. (13) Thus the positive stationary solution (𝑁^∗,𝑃^∗) of system (3) is an unstable equilibrium solution if 𝛿<𝛿_𝐻 and oscillates periodically with respect to time with the frequency 𝜔_𝐻=Im𝜆_𝑘=det(𝐴)=𝛿𝛿−.𝛼+𝛽 (14) The corresponding wavelength is 𝜆_𝐻=2𝜋𝜔_𝐻=2𝜋𝛼+𝛽𝛿.𝛼+𝛽−1 (15)

2.2. Turing Instability

Turing instability (or called Turing bifurcation) is a phenomenon that causes certain reaction-diffusion system to lead to spontaneous stationary configuration. That is why Turing instability is often called diffusion-driven instability. The Turing instability is not dependent on the geometry of the system but only on the reaction rates and diffusion. It can occur only when the inhibitor (𝑃) diffuses faster than the activator (𝑁) [3, 8, 9, 20, 21, 28].

In fact, the Turing instability sets in when at least one of the solutions of (8) crosses the imaginary axis. In other words, the spatially homogeneous steady state will become unstable due to heterogeneous perturbation when at least one solution of (8) is positive. For this reason, at least one out of the following two inequalities is violated to occur the Turing instability phenomenon: tr𝐴_𝑘=𝐹_𝑁+𝐺_𝑃−𝑑₁+𝑑₂𝑘²𝐴<0,det_𝑘=𝑑₁𝑑₂𝑘⁴−𝑑₁𝐺_𝑃+𝑑₂𝐹_𝑁𝑘²+𝐹_𝑁𝐺_𝑃−𝐹_𝑃𝐺_𝑁>0. (16) Since perturbations of zero wavenumber are stable when diffusive instability sets in by definition (due to the stability for the nonspatial steady state), tr(𝐴)=𝐹_𝑁+𝐺_𝑃<0 and det(𝐴)=𝐹_𝑁𝐺_𝑃−𝐹_𝑃𝐺_𝑁>0 are satisfied. It is seen from these facts that the first condition in (16) always holds. Hence we have only one for the instability condition; that is, 𝐻(𝑘²)≡det(𝐴_𝑘)<0. Elementary calculations yield that 𝐻𝑘²=𝑑₁𝑑₂𝑘⁴+𝑑₁𝛿+𝑑₂−𝑑₂−𝑑𝛼+𝛽²𝛽𝛼+𝛽²𝑘²1+𝛿1−𝛼+𝛽<0. (17) Thus the minimum of 𝐻(𝑘²) occurs at the critical wavenumber 𝑘²_𝑇, where 𝑘²_𝑇=𝑑₂𝐹_𝑁+𝑑₁𝐺_𝑃2𝑑₁𝑑₂=12𝑑₁𝑑₂−𝑑₁𝛿−𝑑₂+𝑑₂+𝑑𝛼+𝛽²𝛽𝛼+𝛽². (18) By substituting 𝑘²=𝑘²_𝑇 into 𝐻(𝑘²), we can get a sufficient condition for the Turing instability as follows: 𝐻𝑘²_𝑇=𝐵𝛿²+12𝐷𝐶𝛿+²16𝐵<0, (19) where 𝑑𝐵=−¹4𝑑₂1,𝐶=1−+𝛽𝛼+𝛽𝛼+𝛽²,1𝐷=−1++𝛽𝛼+𝛽𝛼+𝛽². (20) It is from replacing the inequality in (19) by the equality that the critical value of bifurcation parameter 𝛿 for the Turing bifurcation equals 𝛿_𝑇=𝑑₂𝑑₁11−+𝛽𝛼+𝛽𝛼+𝛽²+2𝛽𝛼+𝛽²11−.𝛼+𝛽 (21) At the Turing threshold 𝛿_𝑇, the spatial symmetry of the system is broken and the patterns are stationary in time and oscillatory in space with the wavelength 𝜆_𝑇=2𝜋𝑘_𝑇. (22)

2.3. Wave Bifurcation

The wave instability caused by the wave bifurcation plays an important part in pattern formations in many areas [7]. In fact, mathematically speaking, the wave bifurcation occurs when Im𝜆_𝑘≠0,Re𝜆_𝑘=0at𝑘=𝑘_𝑤≠0. (23) Thus, from elementary calculation, the critical value of wave bifurcation parameter 𝛿 can be obtained as 𝛿_𝑤=𝛽2𝛼+𝛽²+12𝑑²₁𝑑₁+𝑑₂√Φ+2𝑑₁𝑑₂11−+𝛼+𝛽𝛽𝑑²₂𝛼+𝛽², (24) where 𝛽Φ=𝛼+𝛽²𝑑²₁+𝑑²₂𝛽𝛼+𝛽²+2𝑑₁𝑑₂211−−𝛽𝛼+𝛽𝛼+𝛽². (25)

It is well known that, at the wave threshold 𝛿_𝑤, both spatial and temporal symmetries are broken and the patterns are oscillatory in space and time with the wavelength 𝜆_𝑤=2𝜋𝑘_𝑤, (26) where 𝑘²_𝑤=1𝑑₁+𝑑₂1−1−𝛿++𝛽𝛼+𝛽𝛼+𝛽². (27)

If one takes 𝛼=0.4, 𝛽=5/6, and 𝑑₁=1 in system (3), the bifurcation diagram, shown in Figure 1, can be obtained using the above linear stability analysis. Since Turing patterns occur when the inhibitor (𝑃) diffuses faster than the activator (𝑁), we take into account the range of the value 𝑑₂ greater than 𝑑₁=1. According to Figure 1, the Hopf bifurcation line and the Turing bifurcation line intersect at a point (4.0611,0.3587) and the Turing bifurcation line and the wave bifurcation line meet at a point (1,0.0932). In addition, these bifurcation lines split the parametric space into five distinct domains. In domain I, located above two bifurcation lines, the steady state is the only stable solution of system (3). Domain II is the region of the pure Turing instability and domain III is the region of Hopf and wave instabilities. Hopf and Turing instabilities take place in domain IV. On the other hand, all three instabilities occur in domain V. Furthermore, Figure 2 shows the dispersion relations corresponding to several values of two parameters 𝑑₂,𝛿. The transitions of Hopf, Turing, and wave modes from stable to unstable for system (3) can be ascertained clearly via Figure 2. While Hopf is subcritical, Turing and wave bifurcations are supercritical.

Figure 1 

Bifurcation diagrams of system (3) with respect to 𝑑2 and 𝛿 when 𝛼=0.4, 𝛽=5/6, and 𝑑1=1. Hopf bifurcation line: 𝛿=0.3587; Turing bifurcation line: 𝛿=0.0932𝑑2; wave bifurcation line: 𝛿=(7/37)𝑑2+(375/1369)𝑑22−(10√15/1369)375𝑑22−232𝑑2+375+(375/1369). Hopf-Turing bifurcation point: (4.0611,0.3587). Turing-wave bifurcation point: (1,0.0932).

Figure 2 

Dispersion relations: (I) (𝑑2,𝛿)=(2,0.5); (II) (𝑑2,𝛿)=(7,0.5); (III) (𝑑2,𝛿)=(2,0.3); (IV1) (𝑑2,𝛿)=(7,0.2); (IV2) (𝑑2,𝛿)=(2,0.15); (V) (𝑑2,𝛿)=(7,0.3). Solid line: the real parts of the eigenvalues; dashed line: imaginary parts of complex pair of eigenvalues.

3. Spatiotemporal Pattern Analysis via Numerical Simulations

In this section, we will investigate spatiotemporal pattern formations of the spatially extended system (3) in two-dimensional space via numerical examples. In fact, the nonuniform stationary states of system (3) that corresponds to spatial patterns cannot be found analytically. In other words, the analytical methods are not sufficient to fully understand the system, which is a reason why we need to use computer simulations. For this, first, take 𝛼=0.4, 𝛽=5/6, and 𝑑₁=1 in system (3). Then the spatially homogeneous equilibrium (𝑁^∗,𝑃^∗)=(0.1892,0.2270) is always the only nontrivial positive stationary state for the nonspatial system of system (3) regardless of the values 𝑑₂ and 𝛿.

In order to solve partial differential equations numerically, one has to discretize the space and time of the given problem. For this reason, the discrete domain for the continuous domain Ω in system (3) is assumed to be the 200×200 lattice sites and the spacing between the lattice points is assumed to be Δ𝑥=Δ𝑦=1. The time evaluation is also discrete; that is, the time goes in steps of Δ𝑡. And all our numerical simulations employ the zero-flux boundary condition.

System (3) is solved numerically by using a finite difference scheme for the spatial derivatives and an explicit Euler method for the time integration with the time step Δ𝑡=0.01. These numerical schemes for the diffusion equation (3) give stable solutions as long as the equalities 𝑑₁Δ𝑡/(Δ𝑥)²,𝑑₂Δ𝑡/(Δ𝑦)²<1/2 are satisfied (see [31]).

It is well known that the spatiotemporal dynamics of a diffusion-reaction system depends on the choice of initial conditions [3, 9, 32, 33]. In this section, an initial condition is taken as a small amplitude random perturbation around the steady state (𝑁^∗,𝑃^∗) since it is very natural from the biological point of view. We stop the simulation when the numerical solutions either reach a stationary state or show oscillatory behaviors.

In the numerical simulations, different types of spatiotemporal dynamics are observed and we have found that the distributions of predator and prey are always of the same type. Consequently, we can restrict our analysis of pattern formations to one distribution. In this section, we show the distribution of predator, for instance.

Now, we will classify spatiotemporal patter formations of system (3) for each bifurcation domain in Figure 1 by taking advantage of numerical simulations.

As mentioned in Section 2, the domain I in Figure 1 is the region where the steady state is the only stable solution of system (3). Thus we will omit to mention dynamical behaviors of system (3) for this domain.

First, we will investigate dynamical behaviors of system (3) with a random small perturbation initial condition of the stationary solution (𝑁^∗,𝑃^∗) by taking into account the parameters in domain II. As mentioned in Section 2, the domain II is the region of the pure Turing instability. If we take the parameter values 𝑑₂=7 and 𝛿=0.6, we can see from Figure 3 that the random initial distribution (see Figure 3(a)) leads to the formation of a regular macroscopic spotted pattern which prevails over the whole domain at last, and after that the dynamics of the system does not change their behaviors anymore.

(a)

(b)

(c)

(d)

Figure 3 

Snapshots of contour pictures of the time evolution of predator in system (3) when 𝑑2=7 and 𝛿𝐻<𝛿=0.6<𝛿𝑇: (a) 0 iteration, (b) 30000 iterations, (c) 100000 iterations, and (d) 200000 iterations.

From now on, through this section, we will not display the snapshots of the initial pattern as Figure 3(a) since all following figures are obtained with a random small perturbation initial condition of the stationary solution as the previous case.

If we take the parameter values 𝑑₂=7 and 𝛿=0.535, we can observe from Figure 4 that the steady state of spotted pattern and the stripelike pattern coexist. In addition, for the parameters 𝑑₂=7 and 𝛿=0.4, Figure 5 shows that stripelike spatial patterns are prevalent in the whole domain eventually. In fact, we can see that the stripelike patterns grow steadily with time and finally they overwhelm the space domain (see Figure 5(c)) after the spotted pattern and the stripelike pattern coexist (see Figures 5(a) and 5(b)).

(a)

(b)

(c)

Figure 4 

Snapshots of contour pictures of the time evolution of predator in system (3) when 𝑑2=7 and 𝛿𝐻<𝛿=0.535<𝛿𝑇: (a) 30000 iterations, (b) 100000 iterations, and (c) 200000 iterations.

(a)

(b)

(c)

Figure 5 

Snapshots of contour pictures of the time evolution of predator in system (3) when 𝑑2=7 and 𝛿𝐻<𝛿=0.4<𝛿𝑇: (a) 20000 iterations, (b) 40000 iterations, and (c) 200000 iterations.

Next, in order to study dynamical behaviors of system (3) in domain V, set the parameters 𝑑₂=7 and 𝛿=0.34. We can know from Figure 1 that, in this case, all three instabilities occur at the same time. In this context, we display Figure 6 to show snapshots of predator spatial patterns at 50000, 100000, and 200000 iterations. Since Turing instabilities take place in domain V, we can see the stripelike spatial patterns as shown in Figure 5. However, an oscillatory pattern in time can be observed even though we cannot observe this phenomenon via snapshots of Figure 6. To make it clear, we exhibit the local phase portrait of system (3) for a fixed point (0.1892,0.2270) in Figure 7. As one sees, thanks to the Hopf instability, a stable cycle can be observed. Moreover, we can ascertain numerically that the frequency of the periodic oscillations in time 𝜔=0.2685 and the corresponding wavelength 𝜆=23.4. Also, it follows from (14) and (15) that the frequency of the periodic oscillations in time is 𝜔_𝐻=0.2536 and the corresponding wavelength 𝜆_𝐻=24.7738.

(a)

(b)

(c)

Figure 6 

Snapshots of contour pictures of the time evolution of predator in system (3) when 𝑑2=7 and 𝛿𝑤<𝛿=0.34<𝛿𝐻: (a) 50000 iterations, (b) 100000 iterations, and (c) 200000 iterations.

Figure 7 

Dynamical behavior of system (3) for a fixed point (0.1892,0.2270) when 𝑑2=7 and 𝛿=0.34: (a) time series of prey, (b) time series of predator, and (c) the stable limit cycle of system (3).

Thus Figures 3, 4, 5, and 6 show typical Turing patterns of predator in system (3) for the fixed value 𝑑₂=7 when the parameter 𝛿 varies. In fact, as the value 𝛿 decreases, the sequence spotted → spot-stripelike mixtures → stripelike pattern is observed.

Now, consider the parameters 𝑑₂=2 and 𝛿=0.34 in domain III. It is from the linear stability analysis in Section 2 that Hopf and wave instabilities occur. As shown in Figure 8, a special pattern may not be formed such as spotted or stripelike patterns since Turing instabilities do not occur in this domain. However, in this case, an oscillatory wave pattern formation is observed. In order to deal with this phenomenon in more detail, we display the local phase portrait of system (3) for a fixed point (0.1892,0.2270) in Figure 9(a), which shows a stable limit cycle, and the population density of predator in Figure 9(b), which is periodic in time and the space-time plots in Figure 9(c), obtained by piling up the lines 𝑦=100 in each snapshot as time progresses. Thus, we can clearly figure out that Hopf and wave instabilities play an important role in forming a wave and periodic oscillatory pattern in space and time as mentioned in Section 2. In fact, numerical calculations yield that the frequency of the periodic oscillations in time 𝜔=0.2474 and the corresponding wavelength 𝜆=25.4. However, the theoretical values of the frequency and the corresponding wavelength are 𝜔_𝐻=0.2536 and 𝜆_𝐻=24.7738, respectively. However, the theoretical values 𝑘_𝑤=0.0789 and 𝜆_𝑤=79.6768 can be obtained from (26) and (27). From these facts, we can surmise that, in this case, Hopf instabilities have an stronger effect on pattern formation in terms of time than wave instabilities.

(a)

(b)

(c)

Figure 8 

Snapshots of contour pictures of the time evolution of predator in system (3) when 𝑑2=2 and 𝛿𝑇<𝛿=0.34<𝛿𝐻: (a) 1000 iterations, (b) 50000 iterations, and (c) 200000 iterations.

(a)

(b)

(c)

Figure 9 

Dynamical behavior of system (3) when 𝑑2=2 and 𝛿=0.34: (a) phase portrait of system (3) for a fixed point (0.1892,0.2270), (b) time series of predator, and (c) space-time plots.

On the contrary, domain IV is the region of Hopf and Turing instabilities according to Figure 1. To discuss dynamical behaviors of system (3) in this domain, we think about two parameter sets as 𝑑₂=2, 𝛿=0.15 and 𝑑₂=7, 𝛿=0.2. As one sees, Figures 10 and 11 show one of the typical spatial patterns called a spiral pattern instead of patterns shown in Turing domains II and V even if Turing instabilities occur in this domain. In fact, in many recent researches [8, 9, 28, 33], spiral patterns have been reported. However, most of them used specific initial conditions which are not the same as the conditions used in this section to generate a spiral pattern. Thus a diffusive system starting with a small amplitude random perturbation around the steady state could have spiral patterns thanks to the cooperations of Turing and Hopf instabilities.

(a)

(b)

(c)

Figure 10 

Snapshots of contour pictures of the time evolution of predator in system (3) when 𝑑2=2 and 𝛿=0.15<𝛿𝑤: (a) 5000 iterations, (b) 150000 iterations, and (c) 200000 iterations.

(a)

(b)

(c)

Figure 11 

Snapshots of contour pictures of the time evolution of predator in system (3) when 𝑑2=7 and 𝛿=0.2<𝛿𝑤: (a) 5000 iterations, (b) 30000 iterations, and (c) 261000 iterations.

4. Conclusion and Discussion

Various spatiotemporal patterns generated by a semi-ratio-dependent predator-prey system with reaction-diffusion are studied theoretically and numerically. The spatial parametric domain is divided into five regions by three bifurcation lines, Hopf, Turing, and wave bifurcation lines, shown in Figure 1. We show that in domains II and V, the sequence spotted → spot-stripelike mixtures → stripelike pattern is emerged as the growth rate of predator decreases since the Turing instability occurs in domain III, an oscillatory wave pattern emerges due to the Hopf and wave instabilities and in domain IV a typical spiral pattern is investigated because of Hopf and Turing instabilities. Moreover, we calculate numerically the values of the frequency and wavelength of the periodic oscillations in time and compare these values to the theoretical values obtained from (14), (15), (26), and (27).

The authors in [3, 9, 32, 33] have investigated the patterns arising from special initial conditions. They indicated that spatiotemporal dynamics of a reaction-diffusion system depend on the choice of initial conditions. For discussing the evolution process of the spatial pattern formation of system (3) with some specific initial conditions, we employ two specific initial conditions which were considered in [3, 9, 32, 33].

Firstly, we employ an initial condition that predators are introduced in a small, localized region of the circle domain while the population of prey is fixed: 𝑁𝑃=()𝑥,𝑦,0=1.0,𝑥,𝑦,00.2,𝑥−100²+𝑦−100²<100,0,otherwise. (28)(29) Secondly, we take an initial condition that describes prey patch placed into a domain with a constant-gradient predator distribution as follows: 𝑁𝑥,𝑦,0=𝑁^∗−𝜖₁(𝑥−180)(𝑥−720)−𝜖₂,𝑃𝑦−90𝑦−210𝑥,𝑦,0=𝑃^∗−𝜖₃(𝑥−450)−𝜖₄𝑦−135, (30) where 𝜖₁=2×10⁻⁷, 𝜖₂=6×10⁻⁷, 𝜖₃=3×10⁻⁵, and 𝜖₄=6×10⁻⁵.

Figures 12 and 13 illustrate snapshots of contour pictures of the time evolution of the predator in system (3) with the initial conditions (29) and (30), respectively, when 𝑑₂=7 and 𝛿=0.6, for 0,100000, and 1000000 iterations.

(a)

(b)

(c)

Figure 12 

Snapshots of contour pictures of the time evolution of predator in system (3) with the initial condition (29) when 𝑑2=7 and 𝛿=0.6: (a) 0 iteration, (b) 100000 iterations, and (c) 1000000 iterations.

(a)

(b)

(c)

Figure 13 

Snapshots of contour pictures of the time evolution of predator in system (3) with the initial condition (30) when 𝑑2=7 and 𝛿=0.6: (a) 0 iteration, (b) 100000 iterations, and (c) 1000000 iterations.

These figures beg the question that the spotted pattern could be a globally stable pattern even though we do not have any theoretical evidences. Similarly, we can figure out numerically that system (3) with the parameters 𝑑₂=7 and 𝛿=0.3, which lies on domain V, has the same phenomena even if we do not exhibit their snapshots. In this case, the system has a stripelike pattern as shown in Figure 5.

Acknowledgment

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2012-0004725).