#### Abstract

The complex dynamics of two types of tritrophic food chain model systems when the species undergo spatial movements, modeling two real situations of marine ecosystem, are investigated in this study analytically and using numerical simulations. The study has been carried out with the objective to explore and compare the competitive effects of fish and molluscs species being the top predators, when phytoplankton and zooplankton species are undergoing spatial movements in the subsurface water. Reaction diffusion systems have been used to represent temporal evolution and spatial interaction among the species. The two model systems differ in an essential way that the top predators are generalist and specialist, respectively, in two models. “*Wave of Chaos*” mechanism is found to be the responsible factor for the pattern (non-Turing) formation in one dimension seen in the food chain ending with top generalist predator. In the present work we have reported WOC phenomenon, for the first time in the literature, in a three-species spatially extended food chain model system. The numerical simulation leads to spontaneous and interesting pattern formation in two dimensions. Constraints on different parameters under which Turing and non-Turing patterns may be observed are obtained analytically. Diffusion-driven analysis is carried out, and the effect of diffusion on the chaotic dynamics of the model systems is studied. The existence of chaotic attractor and long-term chaotic behavior demonstrate the effect of diffusion on the dynamics of the model systems. It is observed from numerical study that food chain model system with top predator as generalist has very rich dynamics and shows very interesting patterns. An ecosystem having top predator as specialist leads to the stability of the system.

#### 1. Introduction

The interest in the study of chaos in ecological systems has increased in the last couple of decades. Interactions of population in ecological systems are modeled by continuous time models and have been studied extensively in the literature [1]. The species interaction models will be close to reality if the dispersal of the population is also considered, along with the time. Diffusion models form a reasonable basis for studying the interaction of population.

Diffusion is a phenomenon by which the biological population spreads according to the irregular motion of each individual of the population. Diffusion equations have been used effectively to describe the movement of numerous animals in mark-recapture studies [2, 3]. Another aspect of reaction-diffusion theory that has excited numerous applied mathematicians is the realization that simply adding diffusion to certain types of multispecies interactions will cause striking spatial patterns to emerge even in homogeneous environments [4].

In chemistry, seminal work of Turing [5] explained diffusion as a driving force for pattern formation and proposed the idea of Turing instability for the first time in the literature. After his work, many models in chemistry were studied to verify diffusion as the influencing factor for pattern formation, such as the Gierer-Meinhardt model [6], the Sel’kov model [7], the Noyes-Field model used in Belousov-Zhabotinskii reaction [8], and the chemotactic diffusion model [9, 10]. Research in diffusion models has taken into account various factors like stochasticity [11], resource availability, spatial heterogeneity [12, 13], environmental borders [14], the Allee effect [15, 16], evolutionary changes [17], and large-scale phenomena such as weather conditions or dispersal effects [18, 19].

The differential equation models for food chain were investigated by many researchers, and rich dynamical behaviors were found [20–24]. Relatively fewer studies have been done on spatially extended food chain models. Maionchi et al. [25] studied the spatial version of Hasting and Powell model and found that the chaotic attractor of the original model was replaced by a stable fixed point or by a simple attractor. Shen and You [26] performed analytical study of diffusive food chain model to obtain conditions for permanence of the system and extinction of species without any numerical details. Xie [27] studied Turing instability in a diffusive three-species model system.

Two types of aquatic food chain model systems are studied in the presence of diffusion in the present paper. In our earlier work, the nondiffusive forms of the two nonlinear deterministic models of marine ecosystem were studied with a different aim of studying the effect of toxin on population dynamics [28]. However, the realistic model should be a diffusive nonlinear model system in order to consider the movement of the individuals with time. In the first aquatic food chain model system, both specialist and generalist predators are present while the second model has only the specialist predator. As the first step of study, the spatial effect has been incorporated in the form of horizontal (one-dimensional) diffusion in the model systems. The spatial effect has been incorporated, again in both of the model systems in the form of two-dimensional diffusion, as the second step. We explore the emergence of (a) molluscs, zooplankton, and phytoplankton populations and (b) fish, zooplankton, and phytoplankton populations in space and time. We carry out analytical and numerical study of these model systems to observe the effect of diffusion on the dynamics of these model systems.

The objective of this work is to investigate the effects of spatial interaction in marine ecosystem when top predator is either specialist or generalist in two different food chains. The investigation is carried out with the aim to understand the predation mechanism that will stabilize the dynamics of marine environment. This paper is organized as follows. In Section 2, the two diffusive model systems and their origin are discussed. The model systems are analysed for the effect of diffusion, and conditions are derived for emerging patterns in terms of the parameters in Section 3. Section 4 gives the numerical simulation results in both one- and two-dimensional spatial domain. At the end, results are discussed with concluding remarks.

#### 2. Model Systems

The two nonlinear deterministic spatially extended prey-predator model systems studied in this work represent two real situations of marine ecosystem, where the top predator is generalist predator in the first case and the top predator is specialist predator in the second case. In our earlier work [28], we studied the nondiffusive version of the model systems where phytoplankton were toxin producing, and effect of types of toxin release function on the dynamical complexity of the model systems was investigated. The model systems studied in this work are spatially extended and nontoxic version of the model systems studied earlier [28, 29]. The diffusive process which explains the effect of movement should be considered in these model systems because the individual of each population is changing their position with time due to turbulence.

##### 2.1. Model System 1

The model system 1 has both kinds of predators, namely, specialist and generalist. The phytoplankton of density are predated by individuals of zooplankton, the specialist predator, of density , which serves as favorite food for generalist predator molluscs population of density . The interaction of the food-chain model at any location and time satisfies the reaction-diffusion equations given by

where , , , , , , , , , , , , , , and are positive constants. , , and are diffusive constants for phytoplankton, zooplankton, and molluscs, respectively. The various parameters [28, 29] are described here: is the rate of self-growth for phytoplankton, is the rate at which zooplankton will die out in the absence of phytoplankton, ’s are the maximum value that per capita rate can attain, , , , and are the half saturation constants, is the rate of competition among individuals of prey, is the growth rate of molluscs population by sexual reproduction, and , , and are diffusivity constants which signify the spatial movement of phytoplankton, zooplankton, and molluscs populations in ocean due to turbulence.

##### 2.2. Model System 2

The model system 2 describes the interaction of prey, specialist predator, and top specialist predator. At any location in space and time , the phytoplankton population of size is the favorite food of zooplankton of size , which serves as favorite food for the fish population of size , and the food chain is given by reaction diffusion equations:

, , , , , , , , , , , , , , and are positive constants. , , and are diffusive constants for phytoplankton, zooplankton, and fish, respectively. The parameter is the rate of mortality of the specialist predator in the absence of , and is measure of its assimilation efficiency. The other parameters have been explained for model system 1, and the same explanation is valid for this model. The Laplacian operator represents in one-dimensional case in two-dimensional case. The systems are defined on a bounded domain (habitat) and are augmented with appropriate initial conditions and the zero-flux boundary conditions. The diffusive processes present in the marine environment are due to the turbulence in ocean; therefore the values of diffusivity constants are considered numerically the same.

#### 3. Diffusion-Driven Analysis

In this section, we will perform a detailed study of the local dynamics of the two spatial model systems in a two-dimensional domain. The conditions for the Turing instability to occur will be obtained by perturbing the homogeneous steady-state solution. The spatial model systems are studied under positive initial conditions and zero flux boundary conditions, given by where .

##### 3.1. Model System 1

Before proceeding to the diffusion-driven analysis of the model system 1, it is worth discussing in brief the nonspatial model system. The biologically feasible equilibrium points of nonspatial form of model system (1a)–(1c) are , , , and . In , , . In , is a positive root of the quadratic equation , and and are given by and . A detailed discussion and stability analysis of the nonspatial model around these equilibrium points can be seen in the literature [28].

For the linear stability analysis of spatial model (1a)–(1c), it is perturbed with the following two-dimensional spatiotemporal perturbation of the form

where , , are sufficiently small constants, and are the components of wave number along and directions, respectively, and is the wavelength. The system is linearized about the non-trivial interior equilibrium point . The characteristic equation of the linearized version of the spatial model system (1a)–(1c) is given by with where is the wave number given by and is a identity matrix: Our interest is the stability properties of the attracting interior equilibrium point , which will lead to the conditions for Turing instability. From (5) and (6), we get the characteristic equation of the form where ,

The characteristic equation can be rewritten as where An application of the Routh-Hurwitz criteria gives the following theorem.

Theorem 1. *The interior equilibrium point is locally asymptotically stable in the presence of diffusion if and only if the following three conditions are satisfied:*

###### 3.1.1. Diffusive Instability

The spatially homogeneous state will be unstable provided that at least one eigenvalue of characteristic equation (10) is positive. This will occur when at least one of the three inequalities of Theorem 1 does not hold.

For (12a), since , , , and are all positive, the inequality (12a) always holds as from the stability condition of interior point in homogeneous state (cf. [28] for detailed analysis of model in homogeneous state). Therefore, the inequality (12b) can be reversed, that is, where can be rewritten as

The minimum of occurs at given by Consequently, the condition for diffusive instability is .

Theorem 2. *The criterion for diffusive instability for the model system 1 is obtained at the critical wave number of the first perturbations obtained by solving (15).*

##### 3.2. Model System 2

First, we will discuss in brief the nonspatial model system. The biologically feasible equilibrium points of nonspatial form of model system 2 given in (2a)–(2c) are , , , and . In , and . In , is a positive root of the quadratic equations , and . A detailed discussion of nonspatial model can be found in [30, 31].

The linear stability analysis of spatial model (2a)–(2c) is done after perturbing the steady-state with a two-dimensional spatiotemporal perturbation of the form given in (4a)–(4c), used for model system 1. We obtained the same results as in Theorems 1 and 2, with the only difference for model system 2 that takes the form .

#### 4. Numerical Simulations

The numerical simulations are carried out to understand the predation mechanism that will stabilize the dynamics of marine environment and to explore spatiotemporal pattern formation resulting from different types of predation. In this section, we investigate the effect of diffusion on the dynamics of the tritrophic food chain model systems. To explore the spatiotemporal dynamics of the model system, we have numerically solved the system of partial differential equations using finite difference method. Forward difference scheme is used for the reaction part. For diffusion part, central difference scheme is used for one-dimensional case, and standard five-point explicit finite difference scheme is used for two-dimensional diffusion terms. The simulation is carried out at different time levels for both one-dimensional and two-dimensional spatial model systems.

##### 4.1. Model System 1

In one-dimensional case model system 1 takes the form as given below:

The equations (16a)–(16c) are solved numerically with zero flux boundary condition at different time levels , 800, 1000, 5000, and 15000, and results are presented in Figure 1. The parameter values used for Figure 1 are as follows: The system of (16a)–(16c) is numerically solved over 20000 mesh points with spatial resolution and time step . The interior equilibrium for model system 1 is .

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

At time , irregular but less dense dynamics is observed, but as time increases to the irregular chaotic patterns increase. The size of the domain occupied by the irregular chaotic patterns increases and hence occupies the whole region at time , and the chaotic dynamics persists as time increases at . At the same set of parameter values given in (17), the spatial model system (16a)–(16c) displays chaotic attractor shown in Figure 2. We observed that the one-dimensional spatially extended food chain model system 1 displayed irregular chaotic dynamics. Here we observe that the space has alternatively been occupied by regular and chaotic patterns. The size of the domain occupied by the irregular chaotic pattern grows with time, displacing the regular pattern and hence occupying the whole region (cf. Figure 1(b) at location , 8000, 16000 and again check the same locations in Figures 1(c) and 1(d) and finally Figures 1(e) and 1(f)). This phenomenon is called “Wave of Chaos” (WOC) [33, 34].

Figure 2 displays the spatial chaotic attractor. From Figure 2, it is found that the population densities evolve in space leading to chaotic behavior and formation of attractor at . This shows chaotic dynamics of the system in space.

The two-dimensional spatially extended model system 1 takes the form as given below:

The system of (18a)–(18c) is numerically solved over mesh points with spatial resolution and time step . The equilibrium population density for prey, specialist predator and generalist predator is (25.992549, 13.333333, 9.9142). Other parameter values are the same as given in (17). At different time steps, the complex spatial patterns are presented in Figure 3, which were obtained at the parameter values given in (17).

**(a)**

**(b)**

**(c)**

**(d)**

Figure 3(a) describes formation of regular spatial distribution of phytoplankton, zooplankton, and mollusks population densities in ocean at time step, . As time grows, the phytoplankton population density mostly decreases over space, while the zooplankton and mollusks population densities increase resulting in very complex spatial patterns (cf. Figures 3(b), 3(c), and 3(d)). In Figure 3, first, second and third columns, respectively, represent population density distribution of phytoplankton, zooplankton, and mollusks species in marine environment.

The wave of chaos phenomenon, found in one-dimensional space, is responsible for pattern formation seen in two-dimensional spaces. This phenomenon has been observed in two-species spatially extended model systems in the literature by different authors [33, 34]. However for the first time, we have reported WOC phenomenon in a three-species spatially extended food chain model system. Till today, no author has found this phenomenon in any three-species spatially extended model system.

##### 4.2. Model System 2

The model system 2 given in (2a)–(2c) describes the food chain between phytoplankton, zooplankton, and fish population. We study the model system 2 in one dimension. The model system (2a)–(2c) with , and , and takes the following form:

The model system (19a)–(19c) is numerically solved over 20000 mesh points with spatial resolution and . The system is simulated for different time levels , 500, 800, 1000, 5000, and 15000. The parameter values used are The interior equilibrium point is (33.9896, 2.2222, 3.5245). The phytoplankton, zooplankton, and fish population, which are prey, specialist predator, and top specialist predator, respectively, in the model system, displays regular dynamics at (Figure 4(a)). At , onset of irregular dynamics is observed (Figure 4(b)). As time increases to , the population displays chaotic irregular dynamics as seen in Figure 4(c). The chaotic dynamics persists when simulation is done for time, , 5000, and 15000 (Figures 4(d)–4(f)). The spatially extended model (19a)–(19c) displays chaotic attractor given in Figure 5 at the parameter values given in (20). The chaotic attractor has been obtained at time . The spatial evolution of population densities is shown in Figure 5 for model system 2.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

In the food chain ending with fish population as top predator, there has not been found any wave of chaos phenomenon in the one-dimensional spatial domain. The spatial domain shows almost regular patterns over entire space at time and (cf. Figures 4(a) and 4(b)). As we increase the time step for simulation, the population densities in the spatial domain show chaotic patterns at , 1000, 5000, 15000 (cf. Figures 4(c), 4(d), 4(e), and 4(f)). The phytoplankton population here shows a continuous decrease in its density over time. This shows the effect of specialist top predator on the dynamics.

Although, in both of the food chains, one ending with mollusks and another ending with fish, the population density patterns in one dimension are chaotic (cf. Figures 1(f) and 4(f); Figures 2 and 5); however their route to chaos is found to be entirely different.

In two-dimensional case, the model system (2a)–(2c) with , , and and takes the following form:

The system (21a)–(21c) is numerically solved over mesh points with spatial resolution and time step . The equilibrium population density is (33.9896, 2.2222, 3.5245). The system is simulated at time , and the results are presented in Figure 6. Parameter values used in Figure 6 are given in (20). For the spatially extended model system (21a)–(21c), the phytoplankton, zooplankton, and fish population densities patterns are presented in first, second, and third columns of Figure 6, respectively. In this food chain, when the simulation is carried out at different time steps, it is observed that population density distribution evolves into simple spatial patterns (cf. Figures 6(c) and 6(d)). However, simulation at time and gives the initial idea that the population densities will evolve over time into complex spatial patterns.

**(a)**

**(b)**

**(c)**

**(d)**

The spatiotemporal dynamics in the food chain model system which ends with specialist predator is much simpler than the dynamics of food chains ending with generalist predator. This study has investigated and found this interesting result and its possible cause. The numerical simulation shows that Wave of Chaos phenomenon, which is found only in food chain ending with generalist predator, is the possible cause of complex spatiotemporal patterns as seen in model system 1. The absence of WOC phenomenon in model system 2 leads to simple dynamics in spatial distribution of species.

After a detailed analytical and numerical simulation of model systems 1 and 2, in both one and two dimensions, we found that the spatially extended food chain model systems displayed chaotic dynamics leading to interesting spatial patterns. Also, from numerical simulation, we observed that patterns in a food chain with top predator as generalist predator are more complex than the patterns displayed by food chain ending with specialist predator.

#### 5. Discussions and Conclusion

The present study reveals that top specialist predators like fish can stabilize the dynamics of the marine ecosystem. This study suggests that presence of only the generalist predator as top predator in the ecosystem may destabilize the marine environment and can make it chaotic. This supports the study done by many ecologists that impact of more fishing will generally tend to reduce the abundance and diversity of parasites resulting in divergent and complex effects [35].

In the present paper, the pattern formation and spatiotemporal population dynamics in food chain model with diffusion are investigated. From the analytical and numerical results obtained, we found that diffusion plays an important role in the pattern formation of the tritrophic food chain model systems. Complex spatial patterns are obtained due to diffusion of the prey and predator population in a food web model. Further, we obtained chaotic dynamics in the simulation of one-dimensional model system. It is observed that Wave of Chaos (WOC) mechanism is the possible cause for spatiotemporal pattern formation. The spatially extended model systems displayed chaotic attractor. The mathematical conditions for diffusion-driven instability are obtained. Numerical simulation of the two model systems studied in this paper shows that the food chain system with diffusion has rich dynamics including spatial patterns and chaotic dynamics.

We observed that a food chain model system with top predator as generalist as in model system 1 displays more complex spatial patterns (cf. Figure 3). The spatial evolution of the prey and predator population at different time levels displayed chaotic dynamics as the time level is increased (cf. Figures 1(d)–1(f)). However, for model system 2, which is a model for food chain ending with specialist predator, the spatial patterns of prey and predator population density are relatively less complex (cf. Figure 6) and the spatial evolution of the population density is chaotic (cf. Figures 4(d)–4(f)). From the numerical simulation, we found that the dynamical complexity of the tritrophic food chain model ending at generalist predator will pose difficult challenges, if modeled more realistically incorporating the effect of both time and space.

We can conclude from the present study that the presence of fish in marine ecosystem makes the food chain model more stable and helps in making marine environment ecofriendly and nonchaotic whereas effects of overfishing can be divergent and complex. This study supports the ecological study that overfishing can lead to global decline in fish parasites which are critical part of biodiversity, and their loss could substantially impact ecosystem function This study will have a direct impact on fisheries policy [35].

This study reveals that a large variety of different patterns and spatiotemporal dynamics can be found in a food chain system with spatial interaction modeled by self-diffusion. Further research is required to explore the dynamical complexity of spatially extended food chain models with more numbers of species.

#### Acknowledgments

The present research of the author is supported by DST under IU-ATC-Phase 2, Project No. SR/RCUK-DST/IUATC Phase 2/2012-IITM (G) and IIT Mandi under the Project No. IITM/SG/NTK/008. The author would like to thank the referees for their constructive comments, which led to a significant improvement of the original paper.