Next Article in Journal
Further Closed Formulae of Exotic 3F2-Series
Next Article in Special Issue
Between Soft θ-Openness and Soft ω0-Openness
Previous Article in Journal
An Improved Algorithm for Identification of Dominating Vertex Set in Intuitionistic Fuzzy Graphs
Previous Article in Special Issue
Stabilization of the Moving Front Solution of the Reaction-Diffusion-Advection Problem
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Qualitative Analysis of Host–Parasitoid Model with Inclusion of Spatial Refuge Effect

1
Department of Mathematics, University of Poonch Rawalakot, Rawalakot 12350, Pakistan
2
Department of Mathematical Sciences, College of Science, Princess Nourah bint Abdulrahman University, P.O. Box 84428, Riyadh 11671, Saudi Arabia
3
Department of Electricity and Electronics, Institute of Research and Development of Processes, Faculty of Science and Technology, Campus of Leioa (Bizkaia), University of the Basque Country, 48940 Leioa, Spain
*
Author to whom correspondence should be addressed.
Axioms 2023, 12(3), 290; https://doi.org/10.3390/axioms12030290
Submission received: 12 January 2023 / Revised: 2 March 2023 / Accepted: 4 March 2023 / Published: 10 March 2023
(This article belongs to the Special Issue Mathematical Modelling and Applications)

Abstract

:
The objective of this work was to investigate the dynamics of host–parasitoid model with spatial refuge effect. For this, two discrete host–parasitoid models were considered under spatial refuge effect. Suppose that a constant population of hosts may seek refuge and protection from an attack of parasitoids. We found the parametric factors affecting the existence of the equilibrium points and uniqueness of equilibrium points. A local stability analysis of host–parasitoid models was also carried out. Bifurcation theory was used to observe that the host–parasitoid models undergo Neimark–Sacker bifurcation. The effect of the existence of constant refuge effect on the local stability and bifurcation of models was also explored. Hybrid chaos control methodology was used to control the chaotic behavior of model. In addition, numerical simulations, bifurcation diagrams, and phase portraits of the models are also presented.
MSC:
34C60; 37C75; 37G10; 37G35; 39A33; 39A30

1. Introduction

The host–parasitoid model in discrete dynamical systems is a way to describe the population dynamics of a host species and its parasitoid species over discrete-time intervals. This model is based on a set of equations that describe how the populations of the two species change from one time step to the next. The model assumes that the populations of the host and parasitoid species are discrete, meaning they are represented by integers. A model that has secured appreciable recognition from both theoretical and experimental point of view is the host–parasitoid model. The host–parasitoid collaborations are the most accepted topics in mathematical biology as they are necessary to deal with the natural enemy of a pest or insect [1,2,3,4,5,6]. We observe that parasitoids are insect species whose larvae grow as parasites on further insect species. It has been discovered that larvae of parasitoids generally kill and sometimes paralyze its host. However, usually, mature parasitoid individuals are free-living insects. Moreover, parasitoids are characterized, in general, by insects that show one or more larval stages that parasite other arthropods, developing inside them and killing them before the end of their life cycle. Life cycles of parasitoids and their hosts are usually synched, e.g., hosts and parasitoids are univoltine (one generation per year) species. Discrete-time steps corresponding to generations are usually used in such models [7,8,9,10]. Parasitoids commonly include species of wasps, flies, beetles, and worms. On the other hand, a special case is where some preys are fully safe from attack of predators within a progressive or spatial refuge. Spatial refuges can take a variety of forms lying between two extremes, which are given as follows:
(i)
Where the constant ratio of the host populations is safe within refuge.
(ii)
Where the constant number is protected.
The most common among these are the constant proportion refuges. For example, flour moth caterpillars are protected from attack of parasitoids by the ichneumonid (Nemeritis canescen) when they are sufficiently deep in the flour medium to be out of reach of the parasitoid’s ovipositor. In this way, only a proportion of the host’s habitat is searched [11,12,13,14,15]. Before starting our discussion related to the mathematical modeling of host–parasitoid interaction under refuge effects, it is worthwhile to describe the pioneering contribution proposed by Nicholson and Bailey [16]. The mathematical model presented by Nicholson and Bailey is acknowledged as the Nicholson–Bailey (NB) model and is directed by the autonomous planar system of nonlinear difference equations as follows:
H t + 1 = r   H t   e x p ( a   P t )                         P t + 1 = c   H t   ( 1 e x p ( a   P t ) )
where H t and P t are state variables representing population concentrations of host and parasitoid at some instant t , respectively. Moreover, r and c are growth rates of host and parasitoid, respectively, and a denotes the searching efficiency of parasitoids. Further investigation reveals that the NB model has trivial and interior (positive) fixed points. Moreover, a unique positive fixed equilibrium point occurs when r > 1 and is always unstable, which contradicts the fact that the NB model is a noble descriptive of regular host–parasitoid interactions because most of these interactions are stable and do not possess any oscillatory behavior. Various attempts were made by numerous mathematicians and ecologists to formulate appropriate modifications to the classical NB model. One of these modifications is to introduce the refuge effect on the classical NB model. For this, we study two modifications proposed by Hassell in [17]:
Firstly, we consider a constant population of hosts, keeping in view the NB model, to examine the spatial refuge effect. A fraction of host specie α is accessible to the parasitoid specie in every single generation. Consequently, ( 1 α ) is the population within spatial refuge. This is represented by the following model.
H t + 1 = r ( 1 α )   H t + r α H t   e x p ( a   P t ) P t + 1 = α   H t   ( 1 e x p ( a   P t ) )                                    
Relative proportions lead to stability under refuges but large areas of parameter space remain for unstable equilibrium (too few hosts within the refuge) or no equilibrium at all (too many hosts within the refuge).
Secondly, consider constant number refuge, where H 0 hosts are always protected, and we obtain:
H t + 1 = r H 0 + r ( H t H 0 ) e x p ( a   P t ) P t + 1 = ( H t H 0 )   ( 1 e x p ( a   P t ) )    
The spatial refuge can play an important role in the dynamics of the model, potentially leading to the persistence of the host population even in the presence of high parasitoid attack rates.
The main contribution of this article is to explore the dynamics of systems (2) and (3). It would be beneficial first to clarify certain fundamental concepts and mathematical outcomes pertaining to our primary inquiry, prior to delving into an all-encompassing examination of the qualitative behaviors exhibited by these systems. Since systems (2) and (3) are planar discrete-time models governed by autonomous nonlinear systems of difference equations, therefore we start with the dynamical behavior of the general planar systems. The key investigations of the manuscript are synopsized in the following discussion.

2. Stability Analysis around Equilibria

In order to obtain the equilibrium points of the system (2), we transform the system as follows.
x n + 1 = r ( 1 p ) x n + r p x n e a y n , y n + 1 = p x n ( 1 e a y n ) .                                                 ,         0 < p < 1 .
Consider map:
( x y ) ( r ( 1 p ) x + r p x e a y p x ( 1 e a y ) )
Here we introduce,
f ( x , y ) = r ( 1 p ) x + r p x e a y
g ( x , y ) = p x ( 1 e a y )
The equilibria O ( 0 , 0 ) and Ω ( x ,   y ) are obtained, where x = 1 a ln r p 1 r ( 1 p ) r 1 ,   y = 1 a ln r p 1 r ( 1 p ) ,   1 < r < 1 1 p , p ( 0 , 1 ) . To demonstrate the stability (local) analysis of a biologically suitable equilibrium, the following lemma is the best interpretation [18,19].
Lemma 1.
Let  F ( λ ) = λ 2 T r a c e ( J )   λ + D e t e r m i n a n t ( J ) .  This is called the characteristic equation of variational matrix  J  and  λ 1  and    λ 2  are the eigenvalues of this equation. Now, a unique positive equilibrium point is termed as a sink if the absolute values of  λ 1  and   λ 2  is one. Consequently, a sink is always a local asymptotic stable equilibrium point. When  | λ 1 |   >   1  along with  | λ 2 |   > 1 the equilibrium point is known as a repeller; thus, we always have a source which is unstable. When  | λ 1 | < 1  and  | λ 2 | > 1  or  | λ 1 | > 1  and  | λ 2 | < 1 then the equilibrium point is called a saddle and it is always unstable. When the absolute value of either of  λ 1  and    λ 2  is one, then we get a non-hyperbolic equilibrium point.
Taking into account the positivity and existence of model (2) at Ω ( x ,   y ) , one can take r > 1 . In Figure 1, the blue region shows the sink, and the red region shows the source. Moreover, the yellow region shows the Hopf bifurcation. The nonexistence of the model (2) is represented by the white area. Thus, the Jacobian matrix at Ω ( x ,   y ) is:
            J ( x ,   y ) = ( 1 r ( 1 + r p r ) ln ( p r 1 + ( 1 + p ) r ) 1 + r 1 + r r ( 1 + ( 1 + p ) r ) ln ( p r 1 + ( 1 + p ) r ) 1 + r )
and the characteristic polynomial is given as:
F ( λ ) = λ 2 ( 1 + ( 1 + ( 1 + p ) r ) ln ( p r 1 + ( 1 + p ) r ) 1 + r ) λ + r ( 1 + ( 1 + p ) r ) ln ( p r 1 + ( 1 + p ) r ) 1 + r .
By simple calculations, we obtain:
F ( 1 ) = ( 1 + ( 1 + p ) r ) ln [ p r 1 + ( 1 + p ) r ] > 0
F ( 1 ) = 2 + ( 1 + r ) ( 1 + ( 1 + p ) r ) ln [ p r 1 + ( 1 + p ) r ] 1 + r
and
F ( 0 ) = r ( 1 + ( p 1 ) r ) ln ( p r 1 + ( p 1 ) r ) r 1 .
Now, suppose λ 1 , λ 1 be two roots of the system of Equations (2), then the positive equilibrium point is given by Ω ( x ,   y ) = ( 1 a ln r p 1 r ( 1 p ) r 1 ,   1 a ln r p 1 r ( 1 p ) ) .   According to the lemma 1, system (2) has the following conditions related to the stability analysis.
  • The equilibrium point Ω ( x ,   y ) is a sink, if and only if F ( 1 ) > 0 and F ( 0 ) < 1 .
  • The equilibrium point Ω ( x ,   y ) is a source, if and only if F ( 1 ) > 0 and F ( 0 ) > 1 .
  • The equilibrium point Ω ( x ,   y ) is a saddle point, if and only if F ( 1 ) < 0 .
Now, we consider the dynamics of model (3) by replacing H and P with x and y . It is clear that system (3) has two equilibrium points, the trivial equilibrium point (0, 0) and the unique positive equilibrium point ( x ,   y ) = ( x 0 + y 1 e a y , 1 a ln | r ( x x 0 ) x r x 0 |   ) , with r > 1 ,   and r x 0 < x .
The Jacobian matrix computed at ( 0 ,   0 ) is:
J ( 0 ,   0 ) = ( r a r x 0 0 a x 0 ) .
Clearly, ( 0 ,   0 ) is a sink if and only if 0 < r < 1 and a saddle if r > 1 .
The Jacobian matrix at the interior steady state ( x ,   y ) is calculated as:
  J   ( x ,   y ) = ( x r x 0 x x 0 a e a y ( x r x 0 ) y ( e a y 1 ) ( x x 0 ) ( r 1 ) x r ( x x 0 ) a ( 1 + 1 e a y 1 ) ( x r x 0 ) y r ( x x 0 ) )
and the characteristic polynomial is given by:
F ( λ ) = λ 2 ( ( x r x 0 ) ( e a y ( r + a y ) r ) ( e a y 1 ) r ( x x 0 ) ) λ + a e a y ( x r x 0 ) y ( e a y 1 ) ( x x 0 ) .
By simple calculations, we obtain:
{ F ( 1 ) = r ( x 0 + r x 0 2 x ) e a y ( r ( x 0 + r x 0 2 x ) + a ( 1 + r ) ( r x 0 x ) y ) ( e a y 1 ) r ( x x 0 ) F ( 1 ) = ( r 1 ) ( e a y ( r x 0 + a ( x r x 0 ) y ) r x 0 ) ( e a y 1 ) r ( x x 0 ) > 0 F ( 0 ) = a e a y ( x r x 0 ) y ( e a y 1 ) ( x x 0 ) .
From the above mathematical computations, we deduced the following results:
  • The unique steady-state ( x ,   y ) is locally asymptotically stable if and only if:
    a ( 1 + r ) ( r x 0 x ) y r ( x 0 + r x 0 2 x ) < e a y 1   and   a ( x r x 0 ) y ( x x 0 ) < 1 e a y   .
  • The unique steady-state ( x ,   y ) is unstable if and only if:
    1 + a ( 1 + r ) ( r x 0 x ) y r ( x 0 + r x 0 2 x ) < e a y   and   a ( x r x 0 ) y ( x x 0 ) > e a y 1 e a y   .
  • The unique steady-state ( x ,   y ) is saddle if and only if:
    e a y 1 < a ( 1 + r ) ( r x 0 x ) y r ( x 0 + r x 0 2 x ) .
  • The eigenvalues of the characteristic polynomial are complex conjugate with magnitude 1, if and only if:
    | ( x r x 0 ) ( e a y ( r + a y ) r ) ( e a y 1 ) r ( x x 0 ) | < 2 ,
    r = e a y ( x x 0 + e a y ( x 0 + x ( a y 1 ) ) ) a x 0 y .  

3. Neimark–Sacker Bifurcation in Model (2)

Now, we study the presence and trajectory of the Neimark–Sacker bifurcation for the unique positive equilibrium point Ω ( x ,   y ) of model (2). We must find sufficient conditions for the Neimark–Sacker bifurcation. The interior unique positive equilibrium point Ω ( x ,   y ) should be non-hyperbolic for the existence of Neimark–Sacker bifurcation, such that the variational matrix calculated at this point has two complex conjugate eigenvalues with an absolute value equal to one. Discrete-time models corresponding to population can also be investigated it the same way [20,21,22]. The conditions that determine the limitations for the existence of the Neimark-Sacker bifurcation are as follows:
F ( 0 ) = r ( 1 + ( 1 + p ) r ) ln [ p r 1 + ( 1 + p ) r ] 1 + r = 1 ,  
or
ln [ p r 1 + ( 1 + p ) r ] = ( 1 + r ) r ( 1 + ( 1 + p ) r )
and
| ( 1 + ( 1 + ( 1 + p ) r ) ln [ p r 1 + ( 1 + p ) r ] 1 + r ) | < 2 .
Suppose that Ω N S = { ( a , r , p ) Є 3 : r r 0 , | A | < 2 } , then there is change of parameters in a small neighborhood of Ω N S , where Ω N S produces the Neimark–Sacker bifurcation for the positive equilibrium point of system (2). The following two-dimensional map expresses the system of Equation (2) as:
r ˜ ( X Y ) ( r ( 1 p ) X + r p X e a Y p X ( 1 e a Y ) ) .
Let r ˜ indicate the bifurcating parameter, and then the perturbed mapping of the above system is as follows:
( X Y ) ( R ( ( 1 p ) X + p X e a Y ) p X ( 1 e a Y ) )
where R = r 1 + r ˜ and the small perturbation parameter is denoted by | r ˜   | < < 1 . Where x = X 1 a ln r p 1 r ( 1 p ) r 1 and y = Y 1 a ln r p 1 r ( 1 p ) . Then from above map, we have
( X Y ) ( m 11 m 12 m 21 m 22 ) ( x y ) + ( f 1 ( x , y ) f 2 ( x , y ) ) ,
where
f 1 ( x , y ) = m 13 x 2 + m 14 x y + m 15 x 2 y + O ( ( | x | + | y | ) 4 ) ,
f 2 ( x , y ) = m 23 x 2 + m 24 x y + m 25 x 2 y + O ( ( | x | + | y | ) 4 ) ,
m 11 = 1 1 R + ln [ p R 1 + ( 1 + p ) R ] a ;   m 12 = e a ( 1 + ( 1 + ( 1 + e a ) p ) R ) ln [ p R 1 + ( 1 + p ) r ] a ( 1 + r )   ;
m 14 = p ( 1 + e a ( 1 + r p R ) p R ) ( 1 + R ln [ p R 1 + ( 1 + p ) R ] a ( 1 + R ) ) ;   m 21 = 1 1 R + ln [ p R 1 + ( 1 + p ) R ] a ;
m 22 = e a ( 1 + ( 1 + ( 1 + e a ) p ) R ) ln [ p R 1 + ( 1 + p ) R ] a ( 1 + R ) ;   m 24 = p ( 1 + e a ( 1 + R p R ) p R ) ( 1 + R ln [ p R 1 + ( 1 + p ) R ] a ( 1 + R ) ) ;
m 25 = p ( 1 + e a ( 1 + R p R ) p R ) 2 + R ln [ p R 1 + ( 1 + p ) R ] a ( 1 + R ) .
The characteristic equation for the Jacobian matrix at equilibrium point (0, 0) is
λ 2 A ( r ˜ ) λ + B ( r ˜ ) = 0 .
Since ( a , r , p ) Є Ω N S . Hence, the complex conjugate roots can be written as
λ 1 = A ( r ˜ ) i 4 B ( r ˜ ) ( A ( r ˜ ) ) 2 2 ,     λ 2 = A ( r ˜ ) + i 4 B ( r ˜ ) ( A ( r ˜ ) ) 2 2 .
where
A ( r ˜ ) = 1 + ( 1 + ( 1 + p ) R ) ln [ p R 1 + ( 1 + p ) R ] 1 + R
and
B ( r ˜ ) = R ( 1 + ( p 1 ) R ) ln [ p R 1 + ( p 1 ) R ] R 1 .
Thus,
| λ 1 | = | λ 2 | = R ( 1 + ( p 1 ) R ) ln [ p R 1 + ( p 1 ) R ] R 1 .
Now, to obtain the non-degeneracy conditions, we observe that
( d | λ 1 | d r ˜ ) r ˜ = 0 = ( d | λ 2 | d r ˜ ) r ˜ = 0 0 .
Moreover, we have 2 < A ( 0 ) < 2 because ( α , β , r 1 ) Ω N S . Now, consider that
A ( 0 ) = 1 + ( 1 + ( 1 + p ) r 1 ) ln [ p r 1 1 + ( 1 + p ) r 1 ] 1 + r 1 .
Hence, we observe that A ( 0 ) 0 and A ( 0 ) 1 , that is consider A ( 0 ) ± 2 , 0 , 1 .
1 + ( 1 + ( 1 + p ) r 1 ) ln [ p r 1 1 + ( 1 + p ) r 1 ] 1 + r 1 1 1 + ( 1 + ( 1 + p ) r 1 ) ln [ p r 1 1 + ( 1 + p ) r 1 ] 1 + r 1 0 . }
Moreover, B ( 0 ) = r ( 1 + ( 1 + p ) r 1 ) ln [ p r 1 1 + ( 1 + p ) r 1 ] 1 + r 1 .
Make sure that A ( 0 ) ± 2 ,   0 , 1 , and as a result we have λ 1 m ,   λ 2 m 1 , 2 , 3 , 4 when r ˜ = 0 . Henceforth, the roots of the system cannot exist within the intersection of the circle with a unit radius with the synchronized axes if r ˜ = 0 . To obtain the normal form when r ˜ = 0 , we assume that k = A ( 0 ) 2 at ω = 4 B ( 0 ) A ( 0 ) 2 2 . Thus,
k = 1 2 [ 1 + ( 1 + ( 1 + p ) r 1 ) ln ( p r 1 1 + ( 1 + p ) r 1 ) r 1 1 ]
and
ω = 4 r ( 1 + ( 1 + p ) r 1 ) ln [ p r 1 1 + ( 1 + p ) r 1 ] 1 + r 1 ( 1 + ( 1 + ( 1 + p ) r 1 ) ln [ p r 1 1 + ( 1 + p ) r 1 ] 1 + r 1 ) 2 2 .
Furthermore, considering the following transformation:
( x y ) ( m 12 0 k m 11 ω ) ( u ν ) .
Under this transformation, our system can be written as:
( x y ) ( k ω ω k ) ( u ν ) + ( f 1 ˜ ( x , y ) f 2 ˜ ( x , y ) ) .
f 1 ˜ ( x , y ) = m 13 m 12 x 2 + m 14 m 12 x y + m 15 m 12 y 2 + O ( ( | u | + | v | ) 4 )
f 2 ˜ ( x , y ) = ( ( k m 11 ) m 13 m 12 ω m 23 ω ) x 2 + ( ( k m 11 ) m 14 m 12 ω m 24 ω ) x y + ( ( k m 11 ) m 15 m 12 ω m 25 ω ) y 2 + O ( ( | u | + | v | ) 4 ) .
x = m 12 u   and   y = ( k m 11 ) u ω v .
Then consider a real number as follow:
L 1 = ( [ R e ( ( 1 2 λ 1 ) λ 2 2 1 λ 1 τ 20 τ 11 ) 1 2 | τ 11 | 2 | τ 02 | 2 + R e ( λ 2 τ 21 ) ] ) r ˜ = 0 .
where
τ 20 = 1 8 [ f 1 ˜ u u f 1 ˜ v v + 2 f 2 ˜ u v + i ( f 2 ˜ u u f 2 ˜ v v 2 f 1 ˜ u v ) ] ,
τ 11 = 1 4 [ f 1 ˜ u u + f 1 ˜ v v + i ( f 2 ˜ u u f 2 ˜ v v ) ] ,
τ 02 = 1 8 [ f 1 ˜ u u f 1 ˜ v v 2 f 2 ˜ u v + i ( f 2 ˜ u u f 2 ˜ v v + 2 f 1 ˜ u v ) ] ,
τ 21 = 1 16 [ f 1 ˜ u u u + f 1 ˜ u v v + f 2 ˜ u u v + f 2 ˜ v v v + i ( f 2 ˜ u u u + f 2 ˜ u v v f 1 ˜ u u v f 1 ˜ v v v ) ] .
After conducting the aforementioned detailed investigation, we have the following precise result:
Theorem 1.
Assume that (1) holds and  L 1 0 ,  then system (2) undergoes the NeimarkSacker bifurcation at the unique positive equilibrium point  Ω ( x ,   y ) , when the parameter r  varies in a small neighborhood of  r 1 .  Furthermore, if  L 1 < 0 ,  then an attracting invariant closed curve bifurcates from the equilibrium point for  r 1 < r ,  and if  L 1 > 0 , then a repelling invariant closed curve bifurcates from the equilibrium point for  r 1 > r .

4. Neimark–Sacker Bifurcation in Model (3)

The two-dimensional map for model (3) is defined as:
r ˜ ( X Y ) ( R ( X 0 + ( X X 0 ) e a Y ( X X 0 ) ( 1 e a Y ) ) .
Let r ˜ represent the bifurcation parameter where | r ˜ | 1 , and then corresponding perturbed mapping of above system is stated as follows:
( X Y ) ( R ( X 0 + ( X X 0 ) e a Y ) ( X X 0 ) ( 1 e a Y ) ) ,
where R = r 2 + r ˜ and the small perturbation parameter is denoted by | r ˜   | < < 1 . Then afterwards, consider the transformations x = X Υ and y = Y ξ , the above map implies that
( X Y ) ( a 11 a 12 a 21 a 22 ) ( x y ) + ( g 1 ( x , y ) g 2 ( x , y ) ) ,
where
g 1 ( x , y ) = a 13 x 2 + a 14 x y + a 15 x 2 y + O ( ( | x | + | y | ) 4 ) ,
g 2 ( x , y ) = a 23 x 2 + a 24 x y + a 15 x 2 y + O ( ( | x | + | y | ) 4 ) .
a 11 = R ,   a 12 = a e a Y ,   a 21 = 1 e a Y ,   a 22 = a e a Y ( X X 0 ) ,   a 24 = a e a Y .
The characteristic equation for the Jacobian matrix at the equilibrium point (0, 0) is
σ 2 A ( r ˜ ) σ + B ( r ˜ ) = 0 .
Since ( a , r , p ) Ω N S . Hence, the complex conjugate roots can be written as
σ 1 , 2 = A ( r ˜ ) ± i 4 B ( r ˜ ) ( A ( r ˜ ) ) 2 2 .
where
A ( r ˜ ) = e a y ( R + a ( x x 0 ) )
B ( r ˜ ) = a e a y R ( x x 0 ) .  
Thus, | σ 1 | = | σ 2 | = a e a y R ( x x 0 ) . Now, for obtaining the non-degeneracy conditions, we have
( d | σ 1 | d r ) r ˜ = 0 = ( d | σ 2 | d r ) r ˜ = 0 = 4 r 2 a e a y ( x x 0 ) 0 .
Moreover, we have 2 < A ( 0 ) < 2 because ( α , β , r 2 ) Ω N S . Also, A ( 0 ) = a e a y ( x x 0 ) . Now, suppose that A ( 0 ) 0 , 1 that is,
A ( 0 ) = e a y ( r 2 + a ( x x 0 ) ) 0         A ( 0 ) = e a y ( r 2 + a ( x x 0 ) ) 1 . }
To obtain the normal form at r ˜ = 0 , we assume that k = A ( 0 ) 2 at ω = 4 B ( 0 ) A ( 0 ) 2 2 , where
k = e a y ( r 2 + a ( x x 0 ) ) 2     and   ω = 4 a e a y r 2 ( x x 0 ) ( e a y ( r 2 + a ( x x 0 ) ) ) 2 2 .
Furthermore, we consider the following transformation:
( x y ) ( a 12 0 k a 11 ω ) ( u ν ) .
Under this transformation, our system can be written as
( x y ) ( k ω ω k ) ( u ν ) + ( f 3 ˜ ( x , y ) f 4 ˜ ( x , y ) ) .
where
f 3 ˜ ( x , y ) = a 13 a 12 x 2 + a 14 a 12 x y + a 15 a 12 y 2 + O ( ( | u | + | v | ) 4 ) .
f 4 ˜ ( x , y ) = ( a 13 a 12 ω a 23 ω ) x 2 + ( ( k a 11 ) a 14 a 12 ω a 24 ω ) x y + ( ( k a 11 ) a 15 m 12 ω a 25 ω ) y 2 + O ( ( | u | + | v | ) 4 ) .
x = a 12 u and y = ( k a 11 ) u ω v . Now, we consider a real number as follows:
L 2 = ( [ R e ( ( 1 2 σ 1 ) σ 2 2 1 σ 1 ξ 20 ξ 11 ) 1 2 | ξ 11 | 2 | ξ 02 | 2 + R e ( σ 2 ξ 21 ) ] ) r ˜ = 0 .
where
ξ 20 = 1 8 [ f 3 ˜ u u f 3 ˜ v v + 2 f 4 ˜ u v + i ( f 4 ˜ u u f 4 ˜ v v 2 f 3 ˜ u v ) ] ,
ξ 11 = 1 4 [ f 3 ˜ u u + f 3 ˜ v v + i ( f 4 ˜ u u f 4 ˜ v v ) ] ,
ξ 02 = 1 8 [ f 3 ˜ u u f 3 ˜ v v 2 f 4 ˜ u v + i ( f 4 ˜ u u f 4 ˜ v v + 2 f 3 ˜ u v ) ] ,
ξ 21 = 1 16 [ f 3 ˜ u u u + f 3 ˜ u v v + f 4 ˜ u u v + f 4 ˜ v v v + i ( f 4 ˜ u u u + f 4 ˜ u v v f 3 ˜ u u v f 3 ˜ v v v ) ] .
After conducting the detailed investigation described above, we have obtained the following precise result:
Theorem 2.
Assume that map (5) holds and  L 2 0 ,  then system (3) undergoes the NeimarkSacker bifurcation at the unique positive equilibrium point  ( X * ,   Y * ) , when the parameter  r  varies in a small neighborhood of  r 2 .  Furthermore, if  L 2 < 0 ,  then an attracting invariant closed curve bifurcates from the equilibrium point for  r 2 < r ,  and if  L 2 > 0 ,  then a repelling invariant closed curve bifurcates from the equilibrium point for  r 2 > r .

5. Chaos Control

For population models, an important feature considered is controlling chaos and bifurcation. Often, we observe that discrete-time models show more complex behavior than continuous models. Chaos control techniques are used to escape population from unpredictable situations. Here, we use the hybrid control feedback methodology to reduce chaos arising because of bifurcation. We use this technique to reduce chaos emerging because of the Neimark–Sacker bifurcation. As we know that system (2) experiences the Neimark–Sacker bifurcation at the equilibrium point Ω ( x ,   y ) , then a controlled system corresponding to this system can be expressed as [23,24,25,26]:
x n + 1 = α [ r ( 1 p ) x n + r p x n e a y n ] + ( 1 α ) x n , y n + 1 = α [ p x n ( 1 e a y n ) ] + ( 1 α ) y n .                                    
where α ( 0 , 1 ) . Here, the growth rate is increased by α in both populations, the host and parasite, where as ( 1 α ) is a harvesting term which is introduced in a correction to the equations. Moreover, the controlled tactic of system (5) is the fusion of both the feedback control and parameter perturbation. Furthermore, by appropriate selection of the controlled parameter α , the Neimark–Sacker bifurcation for the equilibrium point Ω ( x ,   y ) of the controlled system (5) could be postponed or excluded entirely. Now, the variational matrix of controlled system (5) calculated at equilibrium point Ω ( x ,   y ) is expressed as:
( 1 r ( 1 + r p r ) α ln ( p r 1 + ( 1 + p ) r ) 1 + r ( 1 + r ) α r 1 α + ( 1 + ( 1 + p ) r ) α ln ( p r 1 + ( 1 + p ) r ) 1 + r ) .
Furthermore, the characteristic equation of the Jacobian matrix of the controlled system is
P ( λ ) = λ 2 [ 2 α + ( 1 + ( 1 + p ) r ) α ln ( p r 1 + ( 1 + p ) r ) 1 + r ] λ + 1 α         + ( 1 + ( 1 + p ) r ) α ( 1 + ( 1 + r ) α ) ln ( p r 1 + ( 1 + p ) r ) 1 + r .
The condition for equilibrium point of the controlled system to be locally asymptotically stable is
| 2 α + ( 1 + ( 1 + p ) r ) α ln ( p r 1 + ( 1 + p ) r ) 1 + r |                   < 2 α + ( 1 + ( 1 + p ) r ) α ( 1 + ( 1 + r ) α ) ln ( p r 1 + ( 1 + p ) r ) 1 + r < 2 .
As we know that system (1) experiences the Neimark–Sacker bifurcation at a unique equilibrium point, then a controlled system corresponding to this system can be expressed as:
x n + 1 = α ( r x 0 + r ( x n x 0 ) e x p ( a   y n ) ) + ( 1 α ) x n y n + 1 = α ( ( x n x 0 )   ( 1 e x p ( a   y n ) ) + ( 1 α ) y n        
where α ( 0 ,   1 ) . Moreover, the controlled strategy of the above controlled system is coma bination of both the feedback control and the parameter perturbation. Furthermore, by an appropriate option for the controlled parameter α , the Neimark–Sacker bifurcation for the equilibrium point ( x * , y * ) of controlled system can be postponed or excluded totally. The variational matrix of the controlled system calculated at equilibrium point is:
( 1 + ( 1 + e a y r ) α a e a y r α ( x + x 0 ) α e a y α 1 + α ( 1 + a e a y ( x x 0 ) ) ) .
Additionally, the characteristic equation associated with variational matrix of the controlled system is as follow:
  P ( λ ) = λ 2 ( 2 + α ( 2 + e a y ( r + a x a x 0 ) ) ) λ + ( 1 + α ) 2 + e a y α ( a ( 1 + α ) ( x x 0 ) + r ( 1 + α ( 1 + a x a x 0 ) ) ) .
For the controlled system, the criteria for equilibrium point to be locally asymptotically stable are the following:
| 2 + α ( e a y ( r + a x a x 0 ) 2 ) | < 1 + ( α 1 ) 2 e a y α ( a ( α 1 ) ( x x 0 ) + r ( 1 + α ( a x a x 0 1 ) ) ) < 2 .

6. Numerical Simulations

Here, we discuss the qualitative and complex dynamical behavior of models (2) and (3) with numerical simulations by introducing some examples. In this section, we elaborate our mathematical and theoretical study with the help of numerical simulations by selecting particular parameter values.
Example 1.
Let us consider system (2). Assume that  a = 0.4 , p = 0.85  and  ( x 0 , y 0 ) = ( 8.12 ,   6.36 ) . Then, the system (2) experiences Neimark-Sacker bifurcation when the bifurcation parameter  r 4.614946371656745 .
Furthermore,  ( a , p , r ) = ( 0.4 , 0.85 ,   4.614946371656745 )  system (2) has unique positive steady-state which is represented by  ( x * , y * ) = ( 8.12326 ,   6.3630567 ) . For the parametric values  ( a , p , r ) = ( 0.4 , 0.85 ,   4.614946371656745 )  the characteristic polynomial becomes
P ( λ ) = λ 2 1.2166872417286627 λ + 1 .
Then, the roots of P ( λ ) = 0 , are λ 1 = 0.608343620864331 0.7936737610338864 i and λ 1 = 0.6083436208643314 + 0.7936737610338864 i with | λ 1 | = | λ 2 | = 1 . Therefore, ( a , p , r ) = ( 0.4 , 0.85 ,   4.614946371656745 )   Ω N S , then the corresponding bifurcation diagram for model (2) can be shown in the following figures:
Figure 2 shows that after the bifurcation point r = 4.614946371656745 , the populations described in the model (2) experience quasiperiodic behavior due to emergence of Neimark–Sacker bifurcation. Furthermore, ( 4.614946371656745 ,   6 ] is the non-chaotic region for the model (2), and [ 2 ,   4.614946371656745 ] is its chaotic region. From Figure 2c, it is obvious to observe that the maximum Lyapunov exponents are almost about zero due to the quasiperiodic behavior.
Moreover, it is interesting to see phase portraits of the system (2) for various values of the bifurcation parameter in the chaotic region. For this, we keep fixed a = 0.4 , p = 0.85 and the initial conditions ( x 0 , y 0 ) = ( 8.12 ,   6.36 ) , and choosing some values of r in the chaotic region, then the phase portraits of the system (2) are portrayed in Figure 3 as follows:
Next, we want to observe the effectiveness of the hybrid control method in the chaotic region due to the emergence of the Hopf bifurcation in the model (2). For this, we choose two fixed parameters, that is, a = 0.4 , p = 0.85 , and keep on varying the other two parameters in the approximate intervals. For this, it is suitable to choose the whole chaotic region [ 2 ,   4.614946371656745 ] for the bifurcation parameter r and the whole allowable unit interval [ 0 ,   1 ] for the control parameter α . The stability region (controllability region) of system (6) can be shown by Figure 4, which shows the effectiveness of hybrid chaos control methodology:
Example 2.
Let  a = 0.2 , H 0 = 2.5 , r [ 1.5 ,   6 ]  and  ( x 0 , y 0 ) = ( 16.035 ,   12.402 ) , then system (3) undergoes the NeimarkSacker bifurcation when  r 4.41524 . The bifurcation diagrams for model (3) can be shown in Figure 5:
The Figure 5 shows that at r = 4.41524 and ( x 0 , y 0 ) = ( 16.0356 ,   12.4027 ) , model (3) experiences the Neimark–Sacker bifurcation. Furthermore, at ( a , H 0 , r ) = ( 0.2 ,   2.5 ,   4.41524 )  a unique positive equilibrium point exists. The phase portraits of model (3) for a = 0.2 , H 0 = 2.5 and r [ 1.5 , 6 ] are depicted in Figure 6a–d as:

7. Conclusions

The spatial refuge can have important implications for the population dynamics of host–parasitoid systems. The presence of a refuge can increase the survival of the host population, reduce the parasitoid attack rate, lead to increased competition among parasitoids, and make the system more complex and difficult to predict. Therefore, it is important to consider the effects of spatial structure when modeling or studying host–parasitoid systems.
In this manuscript, two host–parasitoid models are discussed. Both the models are modifications of the Nicholson–Bailey model proposed by Bailey et al., in [16]. One of these modifications is to introduce the refuge effect on the classical NB model. For this, we study the two modifications proposed by Hassell. At first, some qualitative features of host–parasitoid models are discussed. These features include the existence of equilibrium points, the stability of equilibria, and the numerical simulations. Moreover, the mathematical results of [2, 3] have been reviewed. They investigated the outcomes of the existence of refuge on the local stability and bifurcation of the models on the three different host–parasitoid models with the spatial refuge effect. The Taylor series method has been used in proposed model (2). It is also proved that model (3) undergoes the Neimark–Sacker bifurcation around the interior equilibrium point. Moreover, it is observed that both proposed models undergo the Neimark–Sacker bifurcation, and the maximum Lyapunov exponents also show that the models experience the Neimark–Sacker bifurcation. In addition to this, the hybrid chaos control methodology is used to control the chaotic behavior of the proposed models. Numerical simulations are provided for further explanation.

Author Contributions

Conceptualization, M.S.S. and Q.D.; methodology, M.S.S.; software, M.S.S.; validation, M.S.S.; formal analysis, M.S.S. and Q.D.; investigation, M.S.S.; resources, M.S.S.; data curation, M.S.S.; writing—original draft preparation, M.S.S., W.F.A. and M.D.l.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research is funded by the Basque Government under the Grants IT1555-22 and KK-2022/00090; and to MCIN/AEI 269.10.13039/501100011033 for Grant PID2021-358 1235430B-C21; and to MCIN/AEI 269.10.13039/501100011033 for Grant PID2021-358 1235430B-C22. The APC was funded by M. De la Sen.

Data Availability Statement

There is no data associated with the manuscript.

Acknowledgments

Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2023R 371), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia. The author M. De la Sen is grateful to the Basque Government for its support through Grants IT1555-22 and KK-2022/00090; and to MCIN/AEI 269.10.13039/501100011033 for Grant PID2021-358 1235430B-C21; and to MCIN/AEI 269.10.13039/501100011033 for Grant PID2021-358 1235430B-C22.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Agarwal, T. Bifurcation analysis of a host-parasitoid model with the beddington-deangelis functional response. Glob. J. Math. Anal. 2014, 2, 565–569. [Google Scholar] [CrossRef]
  2. Bešo, E.; Mujic, N.; Kalabušić, S.; Pilav, E. Neimark–Sacker bifurcation and stability of a certain class of host-parasitoid models with host refuge effect. Int. J. Bifurc. Chaos 2020, 29, 1950169. [Google Scholar] [CrossRef]
  3. Bešo, E.; Mujic, T.; Kalabušić, S.; Pilav, E. Stability of a certain class of a host-parasitoid model with a spatial refuge effect. J. Biol. Dyn. 2020, 14, 1–31. [Google Scholar] [CrossRef]
  4. Wu, D.; Zhao, H. Global qualitative analysis of a discrete host-parasitoid model with refuge and strong Allee effects. Math. Methods Appl. Sci. 2018, 41, 2039–2062. [Google Scholar] [CrossRef]
  5. Din, Q.; Hussain, M. Controlling chaos and Neimark-Sacker bifurcation in a host-parasitoid model. Asian J. Control. 2019, 21, 1202–1215. [Google Scholar] [CrossRef]
  6. Din, Q.; Saeed, U. Bifurcation analysis and chaos control in a host-parasitoid model. Math. Methods Appl. Sci. 2017, 40, 5391–5406. [Google Scholar] [CrossRef]
  7. Kapçak, S.; Ufuktepe, U.; Elaydi, S. Stability and invariant manifolds of a generalized Beddington host-parasitoid model. J. Biol. Dyn. 2013, 7, 233–253. [Google Scholar] [CrossRef] [PubMed]
  8. Livadiotis, G.; Assas, L.; Dennis, B.; Elaydi, S.; Kwessi, E. A discrete-time host–parasitoid model with at Allee effect. J. Biol. Dyn. 2014, 9, 34–51. [Google Scholar] [CrossRef]
  9. Liu, X.; Chu, Y.; Liu, Y. Bifurcation and chaos in a host-parasitoid model with a lower bound for the host. Adv. Differ. Equ. 2018, 2018, 1–15. [Google Scholar] [CrossRef]
  10. Ma, X.; Din, Q.; Rafaqat, M.; Javaid, N.; Feng, Y. A Density-dependent host-parasitoid Model with stability, bifurcation and chaos control. Mathematics 2020, 8, 536. [Google Scholar] [CrossRef]
  11. Schreiber, S.J. Host-parasitoid dynamics of a generalized Thompson model. J. Math. Biol. 2006, 52, 719–732. [Google Scholar] [CrossRef] [PubMed]
  12. Sebastian, E.; Victor, P. Stability analysis of a host-parasitoid model with logistic growth using Allee effect. Appl. Math. Sci. 2015, 9, 3265–3273. [Google Scholar] [CrossRef]
  13. Wiggins, S. Introduction to Applied Nonlinear Dynamical Systems and Chaos; Springer: New York, NY, USA, 2003. [Google Scholar]
  14. Yuan, L.G.; Yang, Q.G. Bifurcation, invariant curve and hybrid control in a discrete-time predator-prey system. Appl. Math. Model. 2015, 39, 2345–2362. [Google Scholar] [CrossRef]
  15. Luo, X.S.; Chen, G.R.; Wang, B.H.; Fang, J.Q. Hybrid control of period-doubling bifurcation and chaos in discrete nonlinear dynamical systems. Chaos Solitons Fractals 2003, 18, 775–783. [Google Scholar] [CrossRef]
  16. Bailey, V.A.; Nicholson, J. The balance of animal populations. Proc. Zool. Soc. Lond 1935, 105, 551–598. [Google Scholar] [CrossRef]
  17. Hassell, M.P. The Dynamics of Arthropod Predator-Prey Systems; Pritnceton Univeristy Press: Princeton, NJ, USA, 1978. [Google Scholar]
  18. Jury, E.I. Inners and Stability of Dynamic Systems; Research supported by the National Science Foundation; Wiley-Interscience: New York, NY, USA, 1974; Volume 316. [Google Scholar]
  19. Liu, X.; Xiao, D. Complex dynamic behaviors of a discrete-time predator-prey system. Chaos Solitons Fractals 2007, 32, 80–94. [Google Scholar] [CrossRef]
  20. Shabbir, M.S.; Din, Q.; Alabdan, R.; Tassaddiq, A.; Ahmad, K. Dynamical complexity in a class of novel discrete-time predator-prey interaction with cannibalism. IEEE Access 2020, 8, 100226–100240. [Google Scholar] [CrossRef]
  21. Taylor, A.D. Parasitiod competition and the dynamics of host-parasitoid models. Am. Nat. 1988, 132, 417–436. [Google Scholar] [CrossRef]
  22. Tang, S.Y.; Chen, L.S. Chaos in functional response host-parasitoid ecosystem models. Chaos Solitons Fractals 2002, 13, 875–884. [Google Scholar] [CrossRef]
  23. Xu, C.L.; Boyce, M.S. Dynamic complexities in a mutual interference host-parasitoid model. Chaos Solitons Fractals 2005, 24, 175–182. [Google Scholar] [CrossRef]
  24. Jamieson, W.T. On the global behaviour of May’s host-parasitoid model. J. Differ. Equ. Appl. 2019, 25, 583–596. [Google Scholar] [CrossRef]
  25. Tassaddiq, A.; Shabbir, M.S.; Din, Q. Discretization, bifurcation and control for a class of predator-prey interaction. Fractal Fract. 2022, 6, 31. [Google Scholar] [CrossRef]
  26. Kaitala, V.; Ylikarjula, J.; Heino, M. Dynamic complexities in host-parasitoid interaction. J. Theor. Biol. 1999, 197, 331–341. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Topological behavior of equilibrium point where 1 < r < 5 and 0 < p < 1.
Figure 1. Topological behavior of equilibrium point where 1 < r < 5 and 0 < p < 1.
Axioms 12 00290 g001
Figure 2. Bifurcation diagrams and MLE of model (2) with a = 0.4, p = 0.85, and r ∈ [2, 6] (a) bifurcation diagram for host (b) bifurcation diagram for parasitoid (c) MLE.
Figure 2. Bifurcation diagrams and MLE of model (2) with a = 0.4, p = 0.85, and r ∈ [2, 6] (a) bifurcation diagram for host (b) bifurcation diagram for parasitoid (c) MLE.
Axioms 12 00290 g002
Figure 3. Phase portraits of model (2) in xnyn-plane for a = 0.4, p = 0.85 and (x0, y0) = (8.12, 6.36). (a) r = 4.614946371656745 (b) r = 4.6 (c) r = 4.5 (d) r = 3.5.
Figure 3. Phase portraits of model (2) in xnyn-plane for a = 0.4, p = 0.85 and (x0, y0) = (8.12, 6.36). (a) r = 4.614946371656745 (b) r = 4.6 (c) r = 4.5 (d) r = 3.5.
Axioms 12 00290 g003aAxioms 12 00290 g003b
Figure 4. Stability region of the controlled system (5).
Figure 4. Stability region of the controlled system (5).
Axioms 12 00290 g004
Figure 5. Bifurcation diagrams and MLE of model (3) when a = 0.2 and H0 = 2.5. (a) bifurcation for host (b) bifurcation in parasitoid (c) MLE.
Figure 5. Bifurcation diagrams and MLE of model (3) when a = 0.2 and H0 = 2.5. (a) bifurcation for host (b) bifurcation in parasitoid (c) MLE.
Axioms 12 00290 g005
Figure 6. Phase portraits of model (3) in xnyn-plane for a = 0.2, H0 = 2.5. (a) r = 4.41524 (b) r = 4.4 (c) r = 4.3 (d) r = 4.45.
Figure 6. Phase portraits of model (3) in xnyn-plane for a = 0.2, H0 = 2.5. (a) r = 4.41524 (b) r = 4.4 (c) r = 4.3 (d) r = 4.45.
Axioms 12 00290 g006
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Shabbir, M.S.; Din, Q.; Alfwzan, W.F.; De la Sen, M. The Qualitative Analysis of Host–Parasitoid Model with Inclusion of Spatial Refuge Effect. Axioms 2023, 12, 290. https://doi.org/10.3390/axioms12030290

AMA Style

Shabbir MS, Din Q, Alfwzan WF, De la Sen M. The Qualitative Analysis of Host–Parasitoid Model with Inclusion of Spatial Refuge Effect. Axioms. 2023; 12(3):290. https://doi.org/10.3390/axioms12030290

Chicago/Turabian Style

Shabbir, Muhammad Sajjad, Qamar Din, Wafa F. Alfwzan, and Manuel De la Sen. 2023. "The Qualitative Analysis of Host–Parasitoid Model with Inclusion of Spatial Refuge Effect" Axioms 12, no. 3: 290. https://doi.org/10.3390/axioms12030290

APA Style

Shabbir, M. S., Din, Q., Alfwzan, W. F., & De la Sen, M. (2023). The Qualitative Analysis of Host–Parasitoid Model with Inclusion of Spatial Refuge Effect. Axioms, 12(3), 290. https://doi.org/10.3390/axioms12030290

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop