Complex dynamics of a four-species food-web model: An analysis through Beddington-DeAngelis functional response in the presence of additional food

Surbhi Rani; Sunita Gakkhar; Anuraj Singh

doi:10.1515/cmb-2023-0116

Article Open Access

Complex dynamics of a four-species food-web model: An analysis through Beddington-DeAngelis functional response in the presence of additional food

Surbhi Rani , Sunita Gakkhar and Anuraj Singh

Published/Copyright: December 13, 2023

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Computational and Mathematical Biophysics Volume 11 Issue 1

Abstract

The four-dimensional food-web system consisting of two prey species for a generalist middle predator and a top predator is proposed and investigated. The middle predator is predating over both the prey species with a modified Holling type-II functional response. The food-web model is effectively formulated, exhibits bounded behavior, and displays dissipative dynamics. The proposed model’s essential dynamical features are studied regarding local stability. We investigated the four species’ survival and established their persistence criteria. In the proposed model, a transcritical bifurcation occurs at the axial equilibrium point. The numerical simulations reveal the persistence of a chaotic attractor or stable focus. The conclusion is that increasing the food available to the middle predator may make it possible to manage and mitigate the chaos within the food chain.

Keywords: persistence; limit cycles; transcritical bifurcation; chaos

MSC 2010: 34D20; 92B05; 92D25

1 Introduction

Among the most important topics for mathematical biologists and theoretical ecologists is the analysis of the dynamical processes of prey-predator interactions. The prey-predator interaction is an appealing subject of study due to its ubiquitous occurrence and relevance. In recent decades, theoretical ecologists and experimental biologists have focused on prey-predator interactions [13]. Many mathematical models of prey-predator interactions have been developed and investigated to determine various species’ consumption and survival dynamics [10,15]. Ordinary differential equations [14,15], fractional differential equations [5,30], partial differential equations [10,29], stochastic differential equations [16], delay differential equations [32], and difference equations can all be used to model prey-predator interactions. Furthermore, using mathematical models, researchers attempt to detect other environmental impacts, such as the Allee effect, species refuge, harvesting, and pattern development.

Malthus [11] was the first to formulate prey-predator interactions using mathematical modeling in the early nineteenth century. The well-known Lotka-Volterra model was eventually improved by including a logistic growth function for prey species [12], encompassing numerous functional responses and environmental impacts, and these improvements are making prey-predator interactions more realistic [3,9]. In environmental biology and ecology, the dynamics of prey-predator interactions can be significantly influenced by the fear effect [28]. Predator populations have the greatest influence on a prey-predator relationship in direct predation and fear of predation. Many mathematical models have expressed worry about taking direct predation into account.

Because of the possibility of predation, prey species may alter their behavior in the presence of predator species. The mere presence of a predator may have a more significant impact on prey species than direct predation [2,25,31].

Numerical investigations in a tri-trophic food chain model carried out by Hasting and Powell [8] show the existence of chaos. The well-known teacup strange attractor was obtained for a biologically reasonable choice of parameters. The key chaos feature in a nonlinear, coupled, multi-species nonlinear dynamical system is unpredictable behavior. The chaos is sensitively dependent on initial conditions and the choice of parameters. When two predator-prey subsystems in a food chain have oscillations such that their frequencies are not commensurate, the complete system tends to have chaotic dynamics. The food webs of competing species may exhibit chaos in the same way as in food chains. This has been demonstrated in a model of two competing prey and a predator system [6].

Chaos is an interesting dynamic behavior, but its control is a challenge from the resource management point of view, as discussed by Schaffer [22]. The model proposed by Gakkhar and Singh [4] showed the control of chaos when an additional predator is introduced in the usual HP model. The role of additional food for the top predator in a tri-trophic food chain has been investigated by Sahoo and Poria [19]. It was concluded that regular dynamics with additional food can control the chaos. However, he has assumed that the additional food is abundant, and its dynamics are ignored. The impact of additional food on the predator’s survival and persistence in the system is also reported by many investigators [18,20,21,23,24,26,27].

When food availability is limited for the predator, it is observed that the predator may incline toward alternate food (prey) instead of its usual food for survival. The additional food may also be provided for the persistence of predators. These may affect the equilibrium density. Predator-prey models with additional food predict that prey increases the predator population.

In this article, the famous Hastings-Powell food chain [8] has been modified, incorporating additional food for the middle predator population with a Beddington-DeAngelis functional response. Then, investigate the effects of adding additional food to the model. Our main aim is to save the top predator population from extinction in the presence of additional food. We established the feasibility conditions, and local stability of different equilibrium points and numerically explored the chaos control dynamics of the system. The article is organized as follows:

In Section 2, the model system is formulated. In Section 3, positivity, boundedness, and the existence of equilibrium points are discussed. Section 4 is devoted to the local stability of various equilibrium points. In Section 5, the persistence of the system is established. In Section 6, numerical explorations are presented. Finally, we summarize biological indications from our analytical observation, and discussions and conclusions are made in Section 7.

2 Mathematical formulation

Let X and U be the densities of the two prey species. The generalist middle predator, which feeds on both prey species, has a density of Y . Let the top predator, Z , predate on the middle predator, Y . As a result, the following mathematical model for the system’s dynamics is proposed:

(1) d X d T = R 0 X 1 − X K 0 − C 1 A 1 X Y B 1 + B 12 X + B 13 U d U d T = R 1 U 1 − U K 1 − C 2 A 2 U Y B 1 + B 12 X + B 13 U d Y d T = A 1 X Y B 1 + B 12 X + B 13 U + A 2 U Y B 1 + B 12 X + B 13 U − D 1 Y − A 3 Y Z B 2 + Y d Z d T = C 3 A 3 Y Z B 2 + Y − D 2 Z .

The schematic representation of interactions among different species is given in Figure 1.

Figure 1

Schematic representation of interactions among different species.

The positive constants R 0 and R 1 represent the growth rates of prey X and U , respectively, whereas K 0 and K 1 represent carrying capacities. In the absence of food, the constants D i ( i = 1 , 2 ) indicate the loss of predators Y and Z . C i − 1 ( i = 1 , 2 ) represent the conversion rate of prey X and U into predator Y , whereas C 3 represents the conversion rate of Y into Z . B 1 , B 12 , and B 13 are the saturating Beddington-type functional response parameters. B 2 is the middle predator’s half-saturation constant. Since the generalist middle predator Y consumes food from X and U , the Beddington-type functional response is explored for Y . The Holling type-II functional response for top predator Z is considered. When U = 0 , the model (1) changes to the Hastings-Powell food chain. The model (1) includes 16 parameters reduced to ten by inserting the nondimensional variables and parameters listed below:

t = R 0 T , x = X K 0 , u = U K 1 , y = C 1 Y K 0 , z = C 1 Z C 3 K 0 a 1 = A 1 K 0 R 0 B 1 , ω 1 = B 12 K 0 B 1 , ω 2 = B 13 K 0 B 1 α , R = R 1 R 0 , a 2 = C 3 A 3 K 0 C 1 B 2 R 0 , a 3 = C 2 A 2 K 0 C 1 B 1 R 0 , a 4 = A 2 K 0 α B 1 R 0 , d 1 = D 1 R 0 , d 2 = D 2 R 0 , b 2 = K 0 C 1 B 2 .

Accordingly, the nondimensional system takes the form

(2) d x d t = x ( 1 − x ) − a 1 x y 1 + ω 1 x + ω 2 u = F 1 d u d t = R u ( 1 − u ) − a 3 u y 1 + ω 1 x + ω 2 u = F 2 d y d t = a 1 x y 1 + ω 1 x + ω 2 u + a 4 u y 1 + ω 1 x + ω 2 u − d 1 y − a 2 y z 1 + b 2 y = F 3 d z d t = a 2 y z 1 + b 2 y − d 2 z = F 4 .

The initial conditions for the system (2) are as follows:

(3) x ( 0 ) ≥ 0 , y ( 0 ) ≥ 0 , z ( 0 ) ≥ 0 , and u ( 0 ) ≥ 0 .

This article examines the four-species food-web model (2), considering an additional food for the middle predator. The objective is to investigate the effect of additional food on system persistence and chaos control. This model is comprised of two subsystems, both of which are food chains. In the absence of U = 0 or X = 0 , the model reduces to Hasting Powell in the X Y Z and U Y Z subsystems.

3 Preliminaries

This section deals with the positivity and boundedness of the system (2). The positivity ensures the population never goes negative. The boundedness may be interpreted as a natural growth restriction due to limited resources.

3.1 Positive invariance

Let X 1 = x , X 2 = u , X 3 = y , X 4 = z , X = ( X 1 , X 2 , X 3 , X 4 ) T ∈ R 4 , F : R + ⟶ R 4 , and F = ( F 1 , F 2 , F 3 , F 4 ) ∈ C ∞ ( R 4 ) . Then, the system (2) can be written in the matrix form as follows:

X ˙ = F ( X ) = x ( 1 − x ) − a 1 x y 1 + ω 1 x + ω 2 u R u ( 1 − u ) − a 3 u y 1 + ω 1 x + ω 2 u a 1 x y 1 + ω 1 x + ω 2 u + a 4 u y 1 + ω 1 x + ω 2 u − d 1 y − a 2 y z 1 + b 2 y a 2 y z 1 + b 2 y − d 2 z ,

with X ( 0 ) = X 0 .

Since the function is sufficiently smooth and satisfies the Lipschitz condition, the existence and uniqueness theorem guarantees the uniqueness and existence of the initial value problem.

Observe that for X ( 0 ) ∈ R + 4 such that X i = 0 , then X i ( 0 ) ≥ 0 for all ( i = 1 – 4 ) . Accordingly, the solution of the system is positively invariant.

3.2 Boundedness

Theorem 3.1

All the solutions ( x ( t ) , ( y ( t ) , ( z ( t ) ) of the system (2) which initiated in R + 4 are uniformly bounded in the region Γ = ( x , u , y , z ) ∈ R + 4 ; 1 a 4 x ( t ) + 1 a 3 u ( t ) + 1 a 4 y ( t ) + 1 a 4 z ( t ) = μ η + ε , ∀ ε > 0 .

Proof

Define a positive definite function

(4) Ω ( t ) = 1 a 4 x ( t ) + 1 a 3 u ( t ) + 1 a 4 y ( t ) + 1 a 4 z ( t ) .

As Ω ( t ) is differentiable in some maximal interval ( 0 , t b ) for an arbitrary η > 0 , the time derivative of equation (4) along the solution of the system (2) is

d Ω d t + η Ω = x a 4 ( η + ( 1 − x ) ) + u a 3 ( η + R ( 1 − u ) ) + y a 4 ( η − d 1 ) + z a 4 ( η − d 2 ) d Ω d t + η Ω ≤ ( η + 1 ) 2 4 + ( η + R ) 2 4 + ( η − d 1 ) 2 4 + ( η − d 2 ) 2 4 .

Hence, we can find μ > 0 such that

d Ω d t + η Ω ≤ μ ∀ t ∈ ( 0 , t b . )

The theory of differential equation [1,7] gives

0 < Ω ( x , u , y , z ) < μ η ( 1 − e − η t ) + Ω ( x ( 0 ) , u ( 0 ) , y ( 0 ) , z ( 0 ) ) e − η t ∀ t ∈ ( 0 , t b ) ,

and for t b ⟶ ∞ , 0 < Ω ( x , u , y , z ) < μ η . Hence, all the solutions of the initial value problems (2) and (3) remain in R 4 + . The solutions that initiate in R 4 + are confined in the compact region

□ Γ = ( x , u , y , z ) ∈ R + 4 ; 1 a 4 x ( t ) + 1 a 3 u ( t ) + 1 a 4 y ( t ) + 1 a 4 z ( t ) = μ η + ε , ∀ ε > 0 .

4 Existence of equilibrium points

This section briefly summarizes the existence of equilibrium points in the model. The equilibrium points of the model (2) are discussed as follows:

The system (2) has a trivial equilibrium point E 0 = ( 0 , 0 , 0 , 0 ) .
The axial points are E x = ( 1 , 0 , 0 , 0 ) and E u = ( 0 , 1 , 0 , 0 ) .
Remark: The system does not admit equilibrium points on the y and z -axis.
The planar point E x u = ( 1 , 1 , 0 , 0 ) always exists where the two species reach their respective carrying capacities.
The other planar points are E x y = ( x ˜ , 0 , y ˜ , 0 ) , where
x ˜ = d 1 a 1 − d 1 ω 1 , y ˜ = a 1 − d 1 ( 1 + ω 1 ) ( a 1 − d 1 ω 1 ) 2
for
(5) a 1 > d 1 ( 1 + ω 1 ) ,
and E u y = ( 0 , u ˘ , y ˘ , 0 ) , where
u ˘ = d 1 a 4 − d 1 ω 2 , y ˘ = R a 4 ( a 4 − d 1 ( 1 + ω 2 ) ) a 3 ( a 4 − d 1 ω 2 ) 2 ,
with
(6) a 4 > d 1 ( 1 + ω 2 ) .
Remark: No equilibrium point exists on x z , u z , or y z planes.
The unique equilibrium point E 1 = ( x ˆ , 0 , y ˆ , z ˆ ) on u = 0 is obtained by solving the following system of equations:
(7) 1 − x ˆ − a 1 y ˆ 1 + ω 1 x ˆ = 0 a 1 x ˆ 1 + ω 1 x ˆ − d 1 − a 2 z ˆ 1 + b 2 y ˆ = 0 a 2 y ˆ 1 + b 2 y ˆ − d 2 = 0 .
The last equation of equation (7) gives
y ˆ = d 2 ( a 2 − b 2 d 2 ) , a 2 > ( a 1 + b 2 ) d 2 .
We obtain
z ˆ = 1 ( a 2 − b 2 d 2 ) a 1 x ˆ 1 + ω 1 x ˆ − d 1
with the condition x ˆ > x ˜ .

Substitution of y ˆ in the first equation of equation (7) gives the quadratic equation (equation (8)) as follows:

(8) ω 1 ( x ˆ ) 2 + x ˆ ( 1 − ω 1 ) + ( − 1 + a 1 d 2 a 2 − b 2 d 2 ) = 0 .

This gives real, positive, and unique

x ˆ = − ( 1 − ω 1 ) + ( 1 + ω 1 ) 2 − 4 a 1 d 2 ω 1 a 2 − b 2 d 2 2 ω 1 ,

provided a 2 > ( a 1 + b 2 ) d 2 .

As a 2 > ( a 1 + b 2 ) d 2 > d 2 b 2 > 0 , this condition ensures the positivity of y ˆ .

Hence,

x ˆ = − ( 1 − ω 1 ) + ( 1 + ω 1 ) 2 − 4 a 1 d 2 ω 1 a 2 − b 2 d 2 2 ω 1 , y ˆ = d 2 ( a 2 − b 2 d 2 ) , z ˆ = 1 ( a 2 − b 2 d 2 ) a 1 x ˆ 1 + ω 1 x ˆ − d 1 ,

with a 2 > ( a 1 + b 2 ) d 2 , x ˆ > x ˜ .

The another equilibrium point E 2 = ( 0 , u ¯ , y ¯ , z ¯ ) may be obtained from the following system of equations:
(9) R ( 1 − u ¯ ) − a 3 y ¯ 1 + ω 2 u ¯ = 0 a 4 u ¯ 1 + ω 2 u ¯ − d 1 − a 2 z ¯ 1 + b 2 y ¯ = 0 a 2 y ¯ 1 + b 2 y ¯ − d 2 = 0 .
We obtain
u ¯ = − ( 1 − ω 2 ) + ( 1 + ω 2 ) 2 − 4 d 2 a 3 ω 2 R ( a 2 − b 2 d 2 ) 2 ω 2 , y ¯ = d 2 ( a 2 − b 2 d 2 ) , z ¯ = 1 ( a 2 − b 2 d 2 ) a 4 u ¯ 1 + ω 2 u ¯ − d 1
for R a 2 > ( a 3 + R b 2 ) d 2 , u ¯ > u ˘ .
The unique equilibrium point E 3 = ( x ˜ ˜ , u ˜ ˜ , y ˜ ˜ , 0 ) exists when
(10) 0 < a 1 − d 1 ( 1 + ω 1 ) a 1 − d 1 ω 1 < R a 1 a 3 < a 4 − d 1 ω 2 a 4 − d 1 ( 1 + ω 2 ) ,
where
(11) x ˜ ˜ = a 3 ( a 4 − d 1 ω 2 ) − R a 1 ( a 4 − d 1 ( 1 + ω 2 ) ) R a 1 ( a 1 − d 1 ω 1 ) + a 3 ( a 4 − d 1 ω 2 ) , u ˜ ˜ = R a 1 ( a 1 − d 1 ω 1 ) − a 3 ( a 1 − d 1 ( 1 + ω 1 ) ) R a 1 ( a 1 − d 1 ω 1 ) + a 3 ( a 4 − d 1 ω 2 ) , y ˜ ˜ = R ( a 1 + a 4 − d 1 ( 1 + ω 1 + ω 2 ) ) ( a 3 a 4 ( 1 + ω 1 ) + R a 1 2 ( 1 + ω 2 ) − a 1 ( R a 4 ω 1 + a 3 ω 2 ) ) ( R a 1 ( a 1 − d 1 ω 1 ) + a 3 ( a 4 − d 1 ω 2 ) ) 2 .
The point E 4 = ( x * , u * , y * , z * ) denotes unique interior equilibrium point,
x * = 1 − R a 1 ( 1 − u * ) a 3 , y * = d 2 a 2 − b 2 d 2 ,

z * = a 1 1 − R a 1 ( 1 − u * ) a 3 + a 4 u * 1 + ω 1 ( 1 − R a 1 ( 1 − u * ) a 3 ) + ω 2 u * − d 1 1 + b 2 a 2 ,

(12) u * = R C ( 2 R a 1 ω 1 + a 3 ( − 1 − ω 1 + ω 2 ) ) + R a 3 2 C ( R a 2 A 2 − d 2 ( R b 2 A 2 + 4 B ) ) 2 R C B ,
with condition a 2 > b 2 d 2 , R A 2 > 4 B d 2 , u * > 1 − a 3 R a 1 , where A = 1 + ω 1 + ω 2 , B = R a 1 ω 1 + a 3 ω 2 , and C = a 2 − b 2 d 2 .

5 Local stability analysis

At any point ( x , u , y , z ) , the Jacobian matrix of the system (2) is computed as follows:

J = 1 − 2 x − a 1 y + a 1 ω 2 u y ( 1 + ω 1 x + ω 2 u ) 2 a 1 w 2 x y ( 1 + ω 1 x + ω 2 u ) 2 − a 1 x 1 + ω 1 x + ω 2 u 0 a 3 ω 1 u y ( 1 + ω 1 x + ω 2 u ) 2 R − 2 R u − a 3 y + a 3 ω 1 x y ( 1 + ω 1 x + ω 2 u ) 2 − a 3 u 1 + ω 1 x + ω 2 u 0 a 1 y + a 1 ω 2 u y − a 4 ω 1 u y ( 1 + ω 1 x + ω 2 u ) 2 − a 1 ω 2 x y + a 4 y + a 4 ω 1 x y ( 1 + ω 1 x + ω 2 u ) 2 a 1 x + a 4 u 1 + ω 1 x + ω 2 u − d 1 − a 2 z ( 1 + b 2 y ) 2 − a 2 y 1 + b 2 y 0 0 a 2 z ( 1 + b 2 y ) 2 a 2 y 1 + b 2 y − d 2 .

Proposition 5.1

It can be easily observed that the eigenvalues about E o ( 0 , 0 , 0 , 0 ) are 1 , R , − d 1 , − d 2 . Accordingly, E o is a saddle point with an unstable manifold along the x-axis and u-axis.

Theorem 5.2

The equilibrium point E x ( 1 , 0 , 0 , 0 ) is a saddle point.

Proof

The eigenvalues of the Jacobian matrix about E x ( 1 , 0 , 0 , 0 ) are − 1 , R , a 1 − d 1 ( 1 + ω 1 ) 1 + ω 1 , − d 2 .

Hence, E x is a saddle point. The unstable manifold exists along the u -axis. It also has an unstable manifold along y -axis, provided a 1 > d 1 ( 1 + ω 1 ) .

If a 1 < d 1 ( 1 + ω 1 ) , then E x remains a saddle point but has a stable manifold along the y-axis. Furthermore, condition (5) for the existence of E x y is violated in this case.□

Theorem 5.3

The equilibrium point E u ( 0 , 1 , 0 , 0 ) is a saddle point.

Proof

The eigenvalues of the variational matrix about the equilibrium point E u ( 0 , 1 , 0 , 0 ) are 1 , − R , a 4 − d 1 ( 1 + ω 2 ) 1 + ω 2 , − d 2 .

Hence, E u is a saddle point. It has an unstable manifold along the x -axis, and the y -axis provides a 4 > d 1 ( 1 + ω 2 ) .

If a 4 < d 1 ( 1 + ω 2 ) , then E u is also a saddle point but has a stable manifold along the y -axis. Furthermore, it violates the existence condition (6) of E u y .

Hence, E u is saddle under the condition a 4 < d 1 ( 1 + ω 2 ) .

Comparing a 4 < d 1 ( 1 + ω 2 ) with existence condition (6), it is observed that the existence of E u y is possible when E u has an unstable manifold in the y -direction.□

Theorem 5.4

The equilibrium point E x u ( 1 , 1 , 0 , 0 ) is saddle, provided

(13) a 1 − d 1 ( 1 + ω 1 ) + a 4 − d 1 ω 2 > 0 .

Proof

The eigenvalues of the variational matrix about the equilibrium point E x u ( 1 , 1 , 0 , 0 ) are − 1 , − R , a 1 + a 4 1 + ω 1 + ω 2 − d 1 , − d 2 .

Under condition (14), the eigenvalue a 1 + a 4 1 + ω 1 + ω 2 − d 1 is positive.

Hence, E x u is a saddle point with an unstable manifold along the y -axis.

If both conditions (5) and (6) are satisfied and equilibrium points E x y and E u y exist, then E x u is a saddle point. However, the existence of both these equilibrium points is only a sufficient condition for E x u to be a saddle point.□

Theorem 5.5

The equilibrium point E x y ( x ˜ , 0 , y ˜ , 0 ) is locally asymptotically stable provided

(14) ( a 2 − d 2 b 2 ) ( a 1 − d 1 ( 1 + ω 1 ) ) < d 2 ( a 1 − d 1 ω 1 ) 2

(15) R < a 3 a 1 a 1 − d 1 ( 1 + ω 1 ) ( a 1 − d 1 ω 1 )

(16) ( a 1 + d 1 ω 1 ) > ( a 1 − d 1 ω 1 ) ω 1 .

Proof

The Jacobian matrix about the equilibrium point E x y ( x ˜ , 0 , y ˜ , 0 ) gives

a 1 ( a 1 − d 1 ω 1 ) λ 2 + d 1 [ ( a 1 + d 1 ω 1 ) − ( a 1 − d 1 ω 1 ) ω 1 ] λ + d 1 ( a 1 − d 1 ( 1 + ω 1 ) ) ( a 1 − d 1 ω 1 ) = 0 .

The eigenvalues are

λ 1 = a 2 ( a 1 − d 1 ( 1 + ω 1 ) ) ( a 1 − d 1 ω 1 ) 2 + b 2 ( a 1 − d 1 ( 1 + ω 1 ) ) − d 2 , λ 2 = R + a 3 a 1 Λ 1 ,

where Λ 1 = − 1 + d 1 a 1 − d 1 ω 1 .

The eigenvalues λ 1 and λ 2 will be negative for conditions (14) and (15), respectively.

If the equilibrium point E x y exists then the solution of quadratic gives negative eigenvalues provided condition (16) is satisfied.

Hence, E x y is locally asymptotically stable under conditions (14)–(16).

The system will admit periodic solutions in x y -plane when ( a 1 + d 1 ω 1 ) = ( a 1 − d 1 ω 1 ) ω 1 and the equilibrium point becomes nonhyperbolic.□

Theorem 5.6

The equilibrium point E u y ( 0 , u ˘ , y ˘ , 0 ) is locally asymptotically stable, provided

(17) R a 4 ( a 4 − d 1 ( 1 + ω 2 ) ) ( a 2 − d 2 b 2 ) < d 2 a 3 ( a 4 − d 1 ω 2 ) 2 ,

(18) 1 R a 2 < a 4 − d 1 ( 1 + ω 2 ) a 4 − d 1 ω 2 ,

(19) ( a 4 + d 1 ω 2 ) > ( a 4 − d 1 ω 2 ) ω 2 .

Proof

The Jacobian matrix about the equilibrium point E u y gives

a 4 ( a 4 − d 1 ω 2 ) λ 2 + d 1 R [ ( a 4 + d 1 ω 2 ) − ( a 4 − d 1 ω 2 ) ω 2 ] λ + d 1 R ( a 4 − d 1 ( 1 + ω 2 ) ) ( a 4 − d 1 ω 2 ) = 0 .

The values of characteristic roots are

λ 1 = a 2 ( R a 4 ( a 4 − d 1 ( 1 + ω 2 ) ) ) a 3 ( a 4 − d 1 ω 2 ) 2 + b 2 ( R a 4 ( a 4 − d 1 ( 1 + ω 2 ) ) ) − d 2 , λ 1 = 1 + R a 2 Λ 2 ,

where Λ 2 = − 1 + d 1 a 4 − d 1 ω 2 .

Conditions (17) and (18) provide the conditions for the negative value of eigenvalues λ 1 and λ 2 .

If E u y exists then the solution of quadratic gives negative eigenvalues provided condition (19) is satisfied. Hence, E u y is locally asymptotically stable under conditions (17)–(19).

The system will admit periodic solutions in u y -plane when ( a 4 + d 1 ω 2 ) = ( a 4 − d 1 ω 2 ) ω 2 and the equilibrium point becomes nonhyperbolic.□

Remark

The characteristic equation about the equilibrium point E i where i = 1 , 2 , 3 , 4 is given as follows:

λ 4 + A i λ 3 + B i λ 2 + C i λ + D i = 0 .

The Routh-Hurwitz criterion gives the local stability of E i , and it is given as follows:

A i , B i , C i , D i > 0 , A i B i > C i , C i ( A i B i − C i ) − D i A i 2 > 0 .

Since the expressions for the coefficients A i , B i , C i , and D i are complex, they are omitted.

6 Bifurcation

Theorem 6.1

A transcritical bifurcation occurs around the axial equilibrium point E x ( 1 , 0 , 0 , 0 ) in the system (2) if the system parameters satisfy the following condition:

a 1 = a 1 T C = d 1 ( 1 + ω 1 ) ,

where a 1 T C is the critical value of transcritical bifurcation.

Proof

If a 1 = a 1 T C = d 1 ( 1 + ω 1 ) then determinant of Jacobian matrix at E x is zero ( Det ( J E x = 0 ) ). Hence, the Jacobian matrix J E x has a zero eigenvalue.

Let X and Y be the eigenvectors of the matrices J E x at a 1 T C and J E x T at a 1 T C corresponding to zero eigenvalue, respectively, then

X = x 1 x 2 x 3 x 4 = − d 1 0 1 0

and

Y = y 1 y 2 y 3 y 4 = 0 0 1 0 .

Now, using Sotomayor’s theorem [17], we have

Δ 1 = Y T F a 1 ( E x ; a 1 T C ) = 0 0 1 0 0 0 0 0 = 0

and

Δ 2 = Y T [ D F a 1 ( E x ; a 1 T C ) X ] = 0 0 1 0 − 1 0 − − d 1 a 1 0 0 0 0 0 0 0 d 1 a 1 0 0 0 0 0 0 0 0 0 = d 1 a 1 ≠ 0 .

Using the values of x i , y i ( i = 1 – 4 ) , the final expression of Δ 3 is given as follows:

Δ 3 = Y T [ D 2 F a 1 ( E x ; a 1 T C ) ( X , X ) ]

Δ 3 = y 3 [ F 3 x x x 1 2 + F 3 y y x 3 2 + 2 F 3 x y x 1 x 3 ] .

Then,

Δ 3 = − a 1 d 1 ( 1 + ω 1 ) 2 ≠ 0 .

Thus, we observe that Δ 1 = 0 , Δ 2 ≠ 0 , and Δ 3 ≠ 0 .

Hence, according to Sotomayor’s theorem, a transcritical bifurcation occurs in system (2) around the axial equilibrium point E x . The transcritical bifurcation is obtained between the equilibria E x and E x y .□

7 Numerical explorations

Extensive numerical experiments are performed to observe the role of additional food on the dynamics of the tri-trophic food chain. Accordingly, the following set of parametric values are chosen as given by Hasting and Powell [8] for the tri-trophic food chain:

a 1 = 5.0 , ω 1 = 3.0 , a 2 = 0.1 , b 2 = 2.0 , d 1 = 0.4 , d 2 = 0.01 .

The additional parameters are chosen as follows:

R = 1 , a 3 = 0.2 , a 4 = 2.9 , ω 2 = 2.2 .

For this choice of data, Hastings and Powell [8] observed various dynamical behaviors such as stable focus, limit cycle, period-doubling, and chaos with respect to half saturation constant ω 1 . These chaotic dynamics in the phase plane take the shape of a teacup attractor for ω 1 = 3.0 .

The impact of additional food is investigated here by varying the parameter ω 2 , keeping all other parameters fixed. For the gradual increase of the parameter ω 2 , the system (2) switches its stability from chaotic oscillation to limit cycle oscillation, then to stable focus. For ω 2 = 2.2 , the system (2) exhibits chaos (Figure 2). When ω 2 is further increased to ω 2 = 4.9 , the limit-cycle oscillations are observed (Figure 3). The time-series in Figure 4 shows the periodic behavior for ω 2 = 13.8 . However, the stable focus at the interior equilibrium point (0.96483, 0.99859, 0.12500, and 0.42906) is obtained for ω 2 = 13.9 (Figure 5). Thus, the three species Hasting-Powell food chain which was chaotic turned out to have stable dynamics when suitable food is added to the intermediate predator.

$Figure 2 (a) Chaotic solution in x y z xyz phase space and (b) chaotic solution in u y z uyz phase space for ω 2 = 2.2 {\omega }_{2}=2.2 .$

Figure 2

(a) Chaotic solution in x y z phase space and (b) chaotic solution in u y z phase space for ω 2 = 2.2 .

$Figure 3 (a) Periodic limit cycle solution in x y z xyz phase space and (b) periodic limit cycle solution in u y z uyz phase space for ω 2 = 4.9 {\omega }_{2}=4.9 .$

Figure 3

(a) Periodic limit cycle solution in x y z phase space and (b) periodic limit cycle solution in u y z phase space for ω 2 = 4.9 .

$Figure 4 Time series of the system (2) for ω 2 = 13.8 {\omega }_{2}=13.8 .$

Figure 4

Time series of the system (2) for ω 2 = 13.8 .

$Figure 5 Time series of the system (2) for ω 2 = 13.9 {\omega }_{2}=13.9 .$

Figure 5

Time series of the system (2) for ω 2 = 13.9 .

These observations indicate that the additional prey could be used as a biological control parameter for controlling chaotic dynamics.

We consider the following data:

a 1 = 5.0 , ω 1 = 3.0 , a 2 = 0.08 , b 2 = 2.0 , d 1 = 0.4 , d 2 = 0.01 .

In the absence of additional prey ( u = 0 ) , the three species food chain x y z has a limit cycle in x y -plane (see Figure 6) while the z -species goes to extinction as a 2 is decreased in this case. The food for the middle predator is insufficient for its survival.

$Figure 6 Limit cycle in x ‒ y x‒y phase plane for a 2 = 0.08 {a}_{2}=0.08 .$

Figure 6

Limit cycle in x ‒ y phase plane for a 2 = 0.08 .

However, by providing additional food ( u ) to the middle predator, persistence is possible in the four-species food-web system. This is exhibited for the following dataset for additional food:

R = 1 , a 3 = 0.2 , a 4 = 0.08 , ω 2 = 2.2 .

The persistence of all four species in the form of chaotic attractors is possible (Figure 7). The system exhibits various dynamical behaviors.

$Figure 7 Chaotic solution in u ‒ y ‒ z u‒y‒z phase space for a 2 = 0.08 {a}_{2}=0.08 and ω 2 = 2.2 {\omega }_{2}=2.2 .$

Figure 7

Chaotic solution in u ‒ y ‒ z phase space for a 2 = 0.08 and ω 2 = 2.2 .

For further increase in ω 2 = 4 , the z -species survive with the additional food, but the main prey ( x ) goes to extinction (Figure 8). This is because the increase in predator ( y ) due to the availability of additional food ( u ) , increases the predation stress on prey ( x ) , leading to its extinction. Furthermore, it is observed that one of the conditions for persistence is violated in this case.

$Figure 8 Chaotic solution in u ‒ y ‒ z u‒y‒z phase space for a 2 = 0.08 {a}_{2}=0.08 and ω 2 = 4 {\omega }_{2}=4 .$

Figure 8

Chaotic solution in u ‒ y ‒ z phase space for a 2 = 0.08 and ω 2 = 4 .

In Figure 9, a bifurcation diagram is plotted with respect to the parameter a 4 . Here, blue and red colors depict prey density. Cyan and purple colors represent the middle predator and top predator, respectively. The parameters, values are same as mentioned above. The splitting in lines represent the stable equilibria to oscillation in magnitude of population and effect reverse for the merging of lines. This ensures the presence of Hopf bifurcation in system (2) with respect to parameter a 4 .

$Figure 9 Bifurcation diagram of system (2) with respect to parameter a 4 {a}_{4} .$

Figure 9

Bifurcation diagram of system (2) with respect to parameter a 4 .

8 Conclusion

A food-web model comprising four species is developed and investigated. The solution is bounded and positively invariant in a closed and bounded domain, reflecting the well-behaved nature of the model. The local stability has been carried out about its various equilibrium points. The survival of all four species is explored, and conditions are established for their persistence. The occurrence of transcritical bifurcation experienced by the system has been discussed analytically. In the context of our study, the occurrence of a transcritical bifurcation may imply that management or conservation efforts need to carefully consider the parameter values to avoid tipping points that lead to undesirable outcomes. For instance, introducing an invasive species or changes in additional food could trigger such shifts, with cascading effects throughout the food web. The numerical simulations reveal rich dynamics in the system, including stable focus, limit cycle, and chaos. Chaos observed in our model involving an additional food underscores the intricate and uncertain dynamics inherent to ecological systems. This chaotic behavior results from intricate, nonlinear interactions among species and can significantly impact ecological aspects such as biodiversity, stability, and managing natural ecosystems. The appropriate additional food may control the chaos in classical Hasting’s model, giving a stable focus where all four species survive.

Funding information: The first author would like to thank the Council of Scientific and Industrial Research, New Delhi, Government of India (File No. 09/143(0902)/2017 EMR-I), for financial support to carry out her research work and the Department of Mathematics, Indian Institute of Technology Roorkee (IIT Roorkee), for providing stimulating scientific environment and resources. The work of Anuraj Singh is supported by a Core Research Grant, Science Engineering Research Board, Govt. of India (CRG/2021/006380).
Conflict of interest: The authors declare that they have no conflict of interest.
Ethical approval: This research did not require ethical approval.

References

[1] Birkhoff, G, & Rota, JC. (1989). Ordinary Differential Equations, New York, NY: John Wiley. Search in Google Scholar

[2] Creel, S., & Christianson, D. (2008). Relationships between direct predation and risk effects. Trends in Ecology & Evolution, 23(4), 194–201. 10.1016/j.tree.2007.12.004Search in Google Scholar PubMed

[3] Debnath, S., Ghosh, U., & Sarkar, S. (2020). Global dynamics of a tritrophic food chain model subject to the Allee effects in the prey population with sexually reproductive generalized-type top predator. Computational and Mathematical Methods, 2(2), e1079. 10.1002/cmm4.1079Search in Google Scholar

[4] Gakkhar, S., & Singh, A. (2012). Control of chaos due to additional predator in the Hastings-Powell food chain model. Journal of Mathematical Analysis and Applications, 385(1), 423–438. 10.1016/j.jmaa.2011.06.047Search in Google Scholar

[5] Ghosh, U., Pal, S., & Banerjee, M. (2021). Memory effect on Bazykin’s prey-predator model: Stability and bifurcation analysis. Chaos, Solitons & Fractals, 143, 110531. 10.1016/j.chaos.2020.110531Search in Google Scholar

[6] Gilpin, M. E. (1979). Spiral chaos in a predator-prey model. The American Naturalist, 113(2), 306–308. 10.1086/283389Search in Google Scholar

[7] Haque, M., & Venturino, E. (2007). An ecoepidemiological model with disease in predator: the ratio-dependent case. Mathematical methods in the Applied Sciences, 30(14), 1791–1809. 10.1002/mma.869Search in Google Scholar

[8] Hastings, A., & Powell, T. (1991). Chaos in a three-species food chain. Ecology, 72(3), 896–903. 10.2307/1940591Search in Google Scholar

[9] Holling, C. S. (1959). The components of predation as revealed by a study of small-mammal predation of the European pine sawfly1. The Canadian Entomologist, 91(5), 293–320. 10.4039/Ent91293-5Search in Google Scholar

[10] Holmes, E. E., Lewis, M. A., Banks, J., & Veit, R. (1994). Partial differential equations in ecology: spatial interactions and population dynamics. Ecology, 75(1), 17–29. 10.2307/1939378Search in Google Scholar

[11] Kot, M. (2001). Elements of Mathematical Ecology. Cambridge: Cambridge University Press. 10.1017/CBO9780511608520Search in Google Scholar

[12] Lotka, A. (1925). Elements of physical biology. Baltimore: Williams and Wilkins. p. 460. Search in Google Scholar

[13] Meng, X., Liu, R., & Zhang, T. (2014). Adaptive dynamics for a non-autonomous Lotka-Volterra model with size-selective disturbance. Nonlinear Analysis: Real World Applications, 16, 202–213. 10.1016/j.nonrwa.2013.09.019Search in Google Scholar

[14] Morozov, A., Petrovskii, S., & Li, B.-L. (2004). Bifurcations and chaos in a predator-prey system with the Allee effect. Proceedings of the Royal Society of London. Series B: Biological Sciences, 271(1546), 1407–1414. 10.1098/rspb.2004.2733Search in Google Scholar PubMed PubMed Central

[15] Murray, J. D. (2001). Mathematical biology II: spatial models and biomedical applications, vol. 3, New York: Springer. Search in Google Scholar

[16] Perc, M., & Szolnoki, A. (2007). Noise-guided evolution within cyclical interactions. New Journal of Physics, 9(8), 267. 10.1088/1367-2630/9/8/267Search in Google Scholar

[17] Perko, L. (1996). Nonlinear systems: Local theory in differential equations and dynamical systems. Newyork: Springer Science and Business Media. pp. 65–178. 10.1007/978-1-4684-0249-0_2Search in Google Scholar

[18] Prasad, B., Banerjee, M., & Srinivasu, P. (2013). Dynamics of additional food provided predator-prey system with mutually interfering predators. Mathematical Biosciences, 246(1), 176–190. 10.1016/j.mbs.2013.08.013Search in Google Scholar PubMed

[19] Sahoo, B., & Poria, S. (2014). The chaos and control of a food chain model supplying additional food to top-predator. Chaos, Solitons & Fractals, 58, 52–64. 10.1016/j.chaos.2013.11.008Search in Google Scholar

[20] Sahoo, B., & Poria, S. (2015). Effects of additional food in a delayed predator-prey model. Mathematical Biosciences, 261, 62–73. 10.1016/j.mbs.2014.12.002Search in Google Scholar PubMed

[21] Samanta, S., Mandal, A. K., Kundu, K., & Chattopadhyay, J. (2014). Control of disease in prey population by supplying alternative food to predator. Journal of Biological Systems, 22(04), 677–690. 10.1142/S0218339014500272Search in Google Scholar

[22] Schaffer, W. M., & Kot, M. (1986). Differential systems in ecology and epidemiology. Chaos, 8, 158–178. 10.1515/9781400858156.158Search in Google Scholar

[23] Singh, A., Tripathi, D., & Kang, Y. (2023). A modified may Holling tanner model: The role of dynamic alternative resources on species’ survival. Journal of Biological Systems, 32, 2450007. 10.1142/S0218339024500074Search in Google Scholar

[24] Tiwari, V., Tripathi, J. P., Jana, D., Tiwari, S. K., & Upadhyay, R. K. (2020). Exploring complex dynamics of spatial predator-prey system: Role of predator interference and additional food. International Journal of Bifurcation and Chaos, 30(7), 2050102. 10.1142/S0218127420501023Search in Google Scholar

[25] Tripathi, D., & Singh, A. (2023). An eco-epidemiological model with predator switching behavior. Computational and Mathematical Biophysics, 11(1), 20230101. 10.1515/cmb-2023-0101Search in Google Scholar

[26] Tripathi, J. P., Abbas, S., & Thakur, M. (2015). Dynamical analysis of a prey-predator model with Beddington-Deangelis type function response incorporating a prey refuge. Nonlinear Dynamics, 80, 177–196. 10.1007/s11071-014-1859-2Search in Google Scholar

[27] Tripathi, J. P., Tripathi, D., Mandal, S., & Shrimali, M. D. (2023). Cannibalistic enemy-pest model: effect of additional food and harvesting. Journal of Mathematical Biology, 87(4), 58. 10.1007/s00285-023-01991-9Search in Google Scholar PubMed

[28] Wang, X., Zanette, L., & Zou, X. (2016). Modelling the fear effect in predator-prey interactions. Journal of Mathematical Biology, 73(5), 1179–1204. 10.1007/s00285-016-0989-1Search in Google Scholar PubMed

[29] Wang, X., & Zhao, M. (2011). Complex dynamics in a ratio-dependent food-chain model with Beddington-Deangelis functional response. Procedia Environmental Sciences, 10, 135–140. 10.1016/j.proenv.2011.09.024Search in Google Scholar

[30] Yin, C., Chen, Y., & Zhong, S.-M. (2014). Fractional-order sliding mode based extremum seeking control of a class of nonlinear systems. Automatica, 50(12), 3173–3181. 10.1016/j.automatica.2014.10.027Search in Google Scholar

[31] Zanette, L. Y., White, A. F., Allen, M. C., & Clinchy, M. (2011). Perceived predation risk reduces the number of offspring songbirds produce per year. Science, 334(6061), 1398–1401. 10.1126/science.1210908Search in Google Scholar PubMed

[32] Zhang, T., & Zang, H. (2014). Delay-induced Turing instability in reaction-diffusion equations. Physical Review E, 90(5), 052908. 10.1103/PhysRevE.90.052908Search in Google Scholar PubMed

Received: 2023-10-03

Revised: 2023-12-05

Accepted: 2023-12-13

Published Online: 2023-12-13

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/cmb-2023-0116

Keywords for this article

persistence; limit cycles; transcritical bifurcation; chaos

Creative Commons

BY 4.0