Influence of impreciseness in designing tritrophic level complex food chain modeling in interval environment

In this paper, we construct a tritrophic level food chain model considering the model parameters as fuzzy interval numbers. We check the positivity and boundedness of solutions of the model system and find out all the equilibrium points of the model system along with its existence criteria. We perform stability analysis at all equilibrium points of the model system and discuss in the imprecise environment. We also perform meticulous numerical simulations to study the dynamical behavior of the model system in detail. Finally, we incorporate different harvesting scenarios in the model system and deploy maximum sustainable yield (MSY) policies to determine optimum level of harvesting in the imprecise environment without putting any unnecessary extra risk on the species toward its possible extinction.


Mathematical modeling in ecology
Ecology is a branch of biological science with regard to interactions among organisms and their biophysical environment, which includes both biotic and abiotic components. Ecological phenomena are studied at many levels and grow up increasingly as technological and environment impacts have grown very fast. Ecosystem ecologists often focus on flow of energy and recycling of nutrients. In the study of ecology the foremost factor is predator-prey interaction, which configures the energy/biomass flow from one trophic level to higher trophic levels and modulates the corresponding population size. The effects of prey responses to predators on the dynamics of the predator-prey interactions are very significant from the biological point of view, and the predator population has a direct or indirect impact on the prey population. After the pioneering work of Lotka [1] and Volterra [2], many complex predator-prey models have been studied to capture and explain the dynamics of prey-predator systems considering various real-life scenarios. Researchers have postulated many hypotheses to functionally signify the coexistence of © The Author(s) 2020. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. the interacting populations in different ambiances [3][4][5][6][7][8][9][10][11][12]. Researchers have paid considerably high attention mostly on two-species or three-species predator-prey models. Hastings and Powell [13] and Alidoust and Ghahfarokhi [14] have investigated a three-species food chain model for chaotic behavior. Roy and Alam [15] have analyzed autonomous and nonautonomous versions of food chain model incorporating intraspecific competition in top predator. Freedman and Waltman [16] reported the conditions of persistence of a three-species predator-prey model. Moura and Silva [17] have explained the food web and ecological models used in the aquatic ecosystems. Different approaches of mathematical modeling and numerical simulation has been carried out by Kademi et al. [18] for employing the microorganism-based probiotics for removal of the prevalent toxins in the food items. Abrams [19] has shown the relationship between food availability and foraging effort of the species in the ecoenvironment. Nath et al. [20] have shown that refugee and Allee effect in the prey species stabilize the chaos in the tritrophic food chain model. Pyramids, cascades, and synthesis of functionality and stability in the food chain were illustrated by Barbier et al. [21]. Gao and Jiang [22] have explained the stationary distribution of the stochastic food chain chemostatic model and its impact on the environment. Kooi and Poggiale [23] have studied singular perturbation and bifurcation of a bitrophic food chain model. Banerjee and Das [24] have shown the impulsive effect on a tritrophic food chain model with mixed functional responses under seasonal perturbations. Castellonas et al. [25] have shown both Hopf and Bautin bifurcations in a tritrophic food chain model. Huang et al. [26] discussed the dynamical behavior of a food chain model with stage structure and time delays. Matouk et al. [27] have focused the dynamical behavior of the fractional-order Hastings-Powell food chain model and its discretization. Liu and Bai [28] analyzed a stochastic tritrophic food chain model with harvesting. Pal et al. [29] have performed stability and bifurcation analysis of a three-species food chain model with delay. Roy and Karmakar [30] and many others have examined the tritropic food chain model with disease in intermediate predator. Curtsdotter et al. [31] discussed the dynamics of food web based on empirical data. Kar et al. [32] discussed the yield management and resilience in a harvested tritrophic food chain model.

Biomathematical modeling in impreciseness environment
Theories of uncertainty play a very crucial role in analysis of the system behavior of various mathematical models in the field of mathematical biology and ecology as uncertainties are naturally and inherently involved in this field. Mathematically, sometimes uncertainties can be described by random variables resulting from systems of stochastic differential equations [33][34][35]. Note that uncertainties may be both stochastic and nonstochastic in nature. However, construction and analysis of an ecological model with uncertainties in the stochastic sense is again far from reality as fuzziness is involved to fit a suitable probability distributive function. Therefore researchers in last few years have given considerably high attention to construction and analysis of ecological models in fuzzy environment to capture, explain, and implement different ecological phenomena in more realistic and meaningful way. For example, Pal et al. [36] studied a predator-prey model with interval biological parameters; De et al. [37] presented stability analysis of combined project of fish, broiler, and ducks in imprecise environment; Zhang and Zhao [38] studied a diffusive predator-prey system with delays and interval biological parameters to obtain a sustainable optimal harvesting policy; Sharma and Samanta [39] presented a two-species com-

Maximum sustainable yield (MSY) and related research works
In an ecological modeling, the term "maximum sustainable yield (MSY)" plays an important role to explore and manage many ecological systems in the stable manner under fluctuating and unpredictable environmental conditions [46]. In the field of ecology, MSY is the largest yield that can be taken out from the stock of species over a long period without putting any unnecessary extra risk on the species. The concept of MSY aims to maintain the population size at the point of maximum growth rate by harvesting the individuals that would normally be added to the population, allowing the population to continue to be productive for a long period of time. Tropically, MSY is measured at half of the carrying capacity when individuals are able to breed to their maximum rate. At this point, there is a surplus of individuals that can be harvested because growth of the population is at its maximum point due to the large number of reproducing individuals. Above this point, density dependent factors increasingly limit breeding until the population reaches carrying capacity. At this point, there are no surplus individuals to be harvested, and yield drops to zero. The terms B MSY and F MSY are associated with MSY policy and stand for biomass at MSY and fishing/harvesting mortality at MSY, respectively. Recently, the concept of MSY is successfully used in many areas like fisheries management, poultry industry, sustainable food-energy transmission in food chain, and so on. Due to the MSY concept, the net increase in the overall productivity of desirable biotic organism keeping the surrounding environment unaffected leads to an environmentally benign state of productivity [47]. Deploying MSY strategy in the mathematical modeling system, the degree of deterministic outcomes is increased exponentially, and the chances of the stochastic outcomes are minimized proportionally.
Several eminent scientists [48][49][50][51][52][53] [54] first observed the maximum sustainable yield of CPUE (Catch-per-unit-effort) on singlespecies population system. Kar and Matsuda [55] discussed implementation of MSY policy for sustainable management of fishery in presence of strong Allee effect. Legovic et al. [56] have introduced the concept of maximum sustainable total yield (MSTY) and discussed the risk of extinction of multiple species from the ecosystem. The impact of maximum sustainable yield in cooperative model is presented by Legovic and Gecek [57], who observed that species having lower biotic potential and net carrying capacity are going to be extinct as they are harvested in maximum level. Matsuda and Abrams [58] have shown that the appropriate MSY policy (which is precisely maximum sustainable revenue for food webs) in multispecies fisheries system helps to maintain the coexistence of all interacting species. Gosh and Kar [59] illustrated precise MSY policy and species extinction in a prey-predator system involving the Holling-Tanner response function.

Motivation and novelties of the work
In mathematical modeling, researchers make several assumptions to capture the reality in a simplified but reasonably meaningful way. In the field of mathematical biology and ecology, interaction among biotic and abiotic organisms has been signified through several functions, where many parameters are involved, and maximum of these parameters have been considered as constant. However, in reality, it is well perceived that uncertainty and impreciseness toward these parameters cannot be ignored due to aggregation of many human and environmental factors. The values of different biological parameters involved in the model of mathematical biology usually depend on the data collected by experts, where lots of impreciseness are involved. Moreover, many biological parameters change with the environment in which the species live. For example, temperature affects the reproduction rate of virus and bacteria, tidal circulation determines the migration rate of aquatic species, seasonal change leads to the migration time of birds, light intensity controls the rate of photosynthesis of plants, and so on. To explore the mechanism of the uncertain parameters affecting population dynamics, appropriate models are desired taking account of the uncertainties and/or fussiness toward its parameters. Although many works have already been done in this field embedding the uncertainty theory, but still there are lots of scope to develop and extend this area. In this paper, we mathematically analyze a tritrophic level food chain model contemplating the model parameters as dynamic interval numbers and deploy MSY policies in the imprecise environment under various harvesting scenarios like prey harvesting only, harvesting of both prey and top-predator, harvesting of both the middle and top predators, and so on.    The interval [α, β] can be represented by a function I n (x, γ ) = (α n (x)) 1-γ (β n (x)) γ for γ ∈ [0, 1]. This function is called an interval-valued function.

Existence and uniqueness of differential equation with interval-valued function
Proof The differential equation can be written as Then the differential equation becomes (using interval arithmetic operation and property) For a fixed n, let us consider the interval-valued function h n (γ ) = a (1-γ ) n b γ n for γ ∈ [0, 1] for an interval δ n ∈ [a n , b n ]. If h n (γ ) is a continuous and strictly increasing function, then

Food chain model
A food chain may be defined as the interrelationship among different biological species among themselves for sustenance in an ecosystem [60]. A food chain is a linear network of links in a food web starting from producer organisms (such as grass or trees, which use radiation from the Sun to make their food) and ending at apex predator species (like grizzly bears or killer whales, etc.). A food chain displays how each living organism gains its food or energy from each other's. A lower animal feeds on green plants, and a higher organism feeds on smaller, and so on. The following flowchart and diagram depicts the food/energy transfer in a food chain.
Classification of food chain: Grazing food chain: This type of food chain often initiates with green plants. A grazing food chain is noticed in aquatic and grassland ecosystems.
Grassland ecosystem: (i) Detritus food chain: Initial step of detritus food chain is a dead decomposed matter, which also depends on the available solar energy. The dead decomposed matter is organic in nature, which is further spitted into simple nutrients by microorganisms such as fungi and bacteria. Generally, a detritus food chain is seen in forest ecosystem.
The flow diagram of the detritus food chain concept could be expressed as Dead organic matter → Detritivores → predators (ii) Parasitic food chain: In this food chain model, either the producer or the consumer is parasitized and the energy transfers to the lower organisms. The transformation of the energy in this food chain model is not significant in nature. The expression for such parasitic mode of a food chain model can be represented as Trees → Fruit eating birds → Lice and bugs → Bacteria and fungi

Mathematical construction of tritrophic food chain model in an imprecise environment
In population dynamics, the very rudimentary model (known as the Lotka-Volterra model) for interacting species can be described as Here x and y denote the densities of prey and predator populations, respectively. In the absence of predator-prey population grow logistically with intrinsic growth rate p 1 (> 0), q 1 is the capture rate of prey by per predator, q 2 is the energy conversion efficiency of predator, and p 2 is the natural death rate of predator. The model system (3.1) is a very basic prey-predator system, which can be considered as a ditrophic food chain model.
In this paper, we are interested in introducing a food chain model of a tritrophic level consisting of prey, intermediate predator, and top predator (or superpredator) in an imprecise environment. Let x(t), y(t), z(t) be their densities at any time t. We make the following assumptions to formulate the model system: (i) In the absence of predator, the prey population x(t) grows logistically with intrinsic growth rate r, and k is the carrying of the environment. (iv) b 2 and γ are the energy conversion coefficient of predation and natural date rate of the top predator z(t). Based on these assumptions, Ghosh and Kar [53] proposed the following food chain model in the deterministic setup: The deterministic tritrophic food chain model (3.2) in the imprecise environment can be changed so that all the coefficients are interval numbers: Again from the second equation of system (3.4) we have Also, from the last equation of system (3.4) we have Hence all solutions of system (3.4) are nonnegative.
Differentiating this function with respect to time, we have . Applying the theory of differential inequalities, we have Hence all solutions of system (3.4) are bounded in R 3 + .

Equilibrium points of the model system along with its feasibility criteria
The model system (3.4) has four equilibrium points, namely the trivial equilibrium point e 0 (0, 0, 0), axial equilibrium point e 1 (k, 0, 0), planar equilibrium point e 2 (x * p1 , y * p1 , 0), and interior equilibrium point e 3 Note that the trivial equilibrium point e 0 (0, 0, 0) of the model system (3.4) always exists without any restriction on parametric space, whereas the axial equilibrium point e 1 (k, 0, 0) is feasible for k > 0, the planar equilibrium point e 2 (x * p1 , y * p1 , 0) is feasible when a   Proof The variational matrix of the model system (3.4) at the equilibrium point e 1 (k, 0, 0) is

Stability analysis of the model system
The eigenvalues of the matrix are λ 1 = -r Proof The variational matrix of the model system (3.4) at e 2 (x * p1 , y * p1 , 0) can be written as One of the eigenvalues of the variational matrix is λ 1 = 0, and the other eigenvalues are solutions of the quadratic equation λ 2 + λ x * 1 y * 1 = 0. It is to be noted that the planar equilibrium point e 2 (x * p1 , y * p1 , 0) is non-Hyperbolic as one Eigen value of Jacobi matrix of the model system at is zero at this equilibrium point.

Lemma 3.6 The interior equilibrium point e
Proof The variational matrix at e 2 (x * p1 , y * p1 , 0) can be written as The characteristic equation of the variational matrix can be written as By the Routh-Hurwith criterion the equilibrium point e 3 (x * I2 , y * I2 , z * I2 ) of system (3.4) is locally asymptotically stable if

MSY policy in an imprecise environment under various harvesting scenarios
Only prey species (x) is being harvested from the system Suppose we are harvesting only the prey species from the system at a rate ex, where the parameter e denotes the harvesting effort. So the model system (3.4) can be rewritten in the following form: (3.5)

Coexistence equilibrium point
The coexistence equilibrium point of (3.5) is denoted by P 1c (x * , y * , z * ), where It is well understood from these equilibrium points that the investigated prey species and the top predator species decrease if the net effort increases proportionally. The yield function than becomes Only middle predator (y) is being harvested from the system Suppose we are harvesting only the middle predator species from the system at rate ey, where the parameter e denoted the harvesting effort. So the model system (3.4) can be rewritten in the following form: The equilibrium biomass of predator in system (3.6) is . The yield function than become sY (e) = e . Here the equilibrium biomass of predator is independent of the effort, but the yield Y (e) is a linear function of effort e.
Both prey (x) and middle predator (y) are being harvested from the system Suppose we are harvesting both the prey and middle predator species from the system at rates, ex and ey, respectively, where the parameter e denotes the harvesting effort. So the model system (3.4) can be rewritten in the following form: (3.7) Now the coexisting equilibrium of system (3.7) can be written as P 2c (x * 1 , y * 1 , z * 1 ), where .
Note that the prey species x * 1 and the top predator species z * 1 decrease with the increase of the harvesting effort e and prey extinction occurs when the harvesting effort reaches the although the middle predator y * 1 is steadily unaffected. The yield of the model system (3.7) at equilibrium is given as Both middle predator (y) and top predator (z) are being harvested from the system Suppose we are harvesting both the middle predator y and top predator z from the system at rates ey and ez, respectively, where the parameter e denotes the harvesting effort. So the model system (3.4) can be rewritten in the following form: (3.8) The coexisting equilibrium of the model system (3.8) can be written as P 3c (x * 2 , y * 2 , z * 2 ), where At the equilibrium point, we clearly observe that the prey species x * 2 and top predator species z * 2 decrease whereas the middle predator species y * 2 increases when the harvesting effort e varies from its lower level to higher level.
The yield of the model system (3.8) at equilibrium is given as

Numerical simulation and discussion
In this section, we perform meticulous numerical simulations to verify and validate the analytical findings of our model system. We have used mathematical software Matlab (2018) and Matcont [61] to numerically approximate the solution of our model system (3.4). We have used inbuilt code ode45, which is based on an explicit Runge-Kutta ( Table 1). To justify the nature of trivial equilibrium point, we have calculated the eigenvalues of the Jacobi matrix at the trivial equilibrium point for different values of the parameter p as shown in the Table 2. We observed that the trivial equilibrium point e 0 (0, 0, 0) of the imprecise model system (3.4) is always unstable in nature, irrespective of the value of the parameter p, which supports our analytical findings that trivial equilibrium point is always unstable (see Lemma 3.3).
Case-II: Study of the nature of axial equilibrium point e 1 (k, 0, 0) In this case, we simulate the model system (3.4) with the set of values of model parameters as reported in Table 3 and choose the value of parameter p of three different levels (p = 0, 0.6, 1) satisfying the condition a   Table 4 The  Table 3 Case-III: Study of the nature of planar equilibrium point e 2 (x * p1 , y * p1 , 0) In this case, we simulate the model system (3.4) with the set of values of model parameters reported in Table 5 and choose the value of parameter p of three different levels (p = 0, 0.6, 1). We have simulated the model system and calculated the eigenvalues of the Jacobi matrix at this equilibrium point. The eigenvalues are shown in Table 6. We observed that in all cases, one of the eigenvalues is zero, so the planar equilibrium point e 2 (x * p1 , y * p1 , 0) is of nonhyperbolic type. Here the system undergoes a saddle-node bifurcation, where a saddle and a node approach each other, coalesce into a single equilibrium (depicted in Fig. 6) and then disappear. Thus the planar equilibrium point e 2 (x * p1 , y * p1 , 0) is always unstable, which supports our analytical findings that the planar equilibrium point is always unstable (see Lemma 3.5). Figure 6 shows the time series plot of the model system in the    Table 7 and choose the value of parameter p of three different levels (p = 0, 0.6, 1), which satisfy the condition specified in Lemma 3.6. We have simulated the model system and calculated the eigenvalues of the Jacobi matrix at this equilibrium point. The eigenvalues are shown in Table 8. Since the real parts of all eigenvalues of the Jacobi    Table 7 matrix are negative, we can conclude that the interior equilibrium point e 3 (x * I2 , y * I2 , z * I2 ) is stable in nature when the model parameters satisfy the condition specified in Lemma 3.6.
Case-V: Study of MSY policy when only prey species is harvested from the system In this case, we simulate the model system (3.5) with the set of values of model parameters reported in Table 7 and choose the value of parameter p of three different levels (p = 0, 0.6, 1). Here we have varied the harvesting effort e and calculated the values of maximum sustainable yield (MSY) and the corresponding level of maximum effort (e msy ). These values are shown in Table 9.
From the graphs of Fig. 10 we observed that both the prey and top-predator biomass decrease when we gradually increase the harvesting effort e, and the rate of decrease in the population biomass depends on the parameter value p (see the different slopes of de-   creasing blue and pink biomass lines for p = 0, 0.6, 1). We further observed that in all cases a unique MSY exists for a particular optimum level of harvesting effort e, which depends on the parametric interval determined through the model parameter p. It is interesting to observe that in case of only prey harvesting from the system, both the MSY and optimum level of harvesting effort increase when the parameter p increases in the interval [0, 1]. Case-VI: Study of MSY policy when both the prey species (x) and middle predator (y) are harvested from the system Here we simulate the model system (3.7) with the set of values of model parameters reported in Table 7 and choose the value of parameter p of three different levels (p = 0, 0.6, 1). We have varied the harvesting effort e in the range [0, 2.5] and calculated the value of maximum sustainable yield (MSY) and the corresponding level of maximum effort (e msy ). These values are shown in Table 10.
From the graphs of Fig. 11 we observed that both the prey and top-predator biomass decrease when we gradually increase the harvesting effort e in the range [0, 2.5] and the rate of decrease in the population biomass depends on the parameter value p (see the different slopes of decreasing blue and pink biomass lines for p = 0, 0.6, 1). We further observed that in all cases a unique MSY exists for a particular optimum level of harvesting effort e, which depends on the parametric interval determined through the model parameter p. It is interesting to observe that in case of harvesting of both the prey species and middle predator from the system, both the MSY and optimum level of harvesting effort increase when the parameter p increases in the interval [0, 1].

Figure 11
The changes in yield and population biomass with respect to the change in harvesting effort e in the range [0, 2.5]. Here the value of parameter p has been varied in three different levels (the leftmost figure is for p = 0; the middle figure is for p = 0.6, and the rightmost figure is for p = 1)

Case-VII: Study of MSY policy when both the middle and top predators are harvested from the system
Here we simulate the model system (3.8) with the set of values of model parameters reported in Table 11 and choose the value of parameter p of three different levels (p = 0, 0.6, 1). We have varied the harvesting effort e in the range [0, 5] and calculated the value of maximum sustainable yield (MSY) and the corresponding level of maximum effort (e msy ). These values are shown in Table 12.
From the graphs of Fig. 12 we observed that both the prey and top-predator biomass decrease whereas middle-predator biomass increases when we gradually increase the harvesting effort e in the range [0,5]. Note that the rate of decrease in the prey population biomass greatly depends on the parameter value p (see the different slopes of decreasing lines of first graph of Fig. 12), whereas the rate of change in predator species (both middle and top predators) with respect to the change of the parameter value p is not much significant (the slopes of different lines of the second and third graphs of Fig. 12 are almost the same). From the fourth graph of Fig. 12 we observed that unique MSY exists for a particular optimum level of harvesting effort e, which depends on the parametric interval determined through the model parameter p.

Figure 12
The changes in the population biomass and total yield with respect to the change in harvesting effort e in the range [0,5]. Here the value of parameter p has been varied in three different levels (p = 0, 0.6, 1)

Conclusion
In this paper, we analyzed a tritrophic level food chain model considering the model parameters as dynamic interval numbers. At first, we have constructed a mathematical model in step-by-step fashion and checked its well-posedness through the positivity and boundedness of solution of the model system. We have found out all the equilibrium points of the model system along with its existence criteria. Stability analysis at all the equilibrium points of the model system was performed, and dynamics of the system was discussed in detail. We have performed meticulous numerical simulations to verify and validate the analytical findings of our model system. Finally, we have incorporated different harvesting scenarios like prey harvesting only, harvesting of both prey and top-predator, harvesting of both the middle and top predators, and so on in the model system and deployed MSY policies in the dynamic interval environment. The model outcome helps to determine the optimum level of harvesting in the imprecise environment without putting any unnecessary extra risk on the species toward its extinction. The study enables the experimentalists to further execute several comprehensive measures while designing inter-species relationbased activities in nature, which would pave the way for securing the endangered and vulnerability of several species of importance, preserving thereby the natural ecosystem as a whole.

Funding
This research was funded by Department of Law, Economics and Human Sciences-Decisions_LAB grant number 2/2020 and the APC was funded by Piano Dipartimento di Eccellenza L.232/2016.

Availability of data and materials
This is not applicable in this paper.