The Spatiotemporal Dynamics of Insect Predator–Prey System Incorporating Refuge Effect

The insect predator–prey system mediates several feedback mechanisms which regulate species abundance and spatial distribution. However, the spatiotemporal dynamics of such discrete systems with the refuge effect remain elusive. In this study, we analyzed a discrete Holling type II model incorporating the refuge effect using theoretical calculations and numerical simulations, and selected moths with high and low growth rates as two exemplifications. The result indicates that only the flip bifurcation opens the routes to chaos, and the system undergoes four spatiotemporally behavioral patterns (from the frozen random pattern to the defect chaotic diffusion pattern, then the competition intermittency pattern, and finally to the fully developed turbulence pattern). Furthermore, as the refuge effect increases, moths with relatively slower growth rates tend to maintain stability at relatively low densities, whereas moths with relatively faster growth rates can induce chaos and unpredictability on the population. According to the theoretical guidance of this study, the refuge effect can be adjusted to control pest populations effectively, which provides a new theoretical perspective and is a feasible tool for protecting crops.

A series of experiments have shown that predator-prey interactions among insects (e.g., moths and spiders) conform to the Holling type II functional response [24][25][26].In most prior studies, the Holling type II predator-prey models were continuous and analyzed by phase plane and bifurcation diagrams, without integrating ecological phenomenon to analyze dynamic behaviors [25][26][27].Although some scholars discretized the Holling type II model, the chaotic phase plane diagrams failed to be further distinguished [28][29][30].These Entropy 2024, 26,196 2 of 21 gaps have largely hindered our explorations of the complexity in insect predator-prey systems, thus further studies on the discrete Holling type II models are urgently needed.
In recent years, extensive studies have incorporated the refuge effect into the predatorprey model, confirming the important role of refuge [27,[31][32][33].The continuous Lotka-Volterra predator-prey models show that the refuge effect has a profound effect on predatorprey interaction [34].Moreover, the predator-prey models with ratio dependence also demonstrate that the refuge effect supports the coexistence of the predator and prey [35].Notably, the discrete-time predator-prey system, which is more compatible in describing the dynamics of insects than the continuous system considering the generation nonoverlapping of many insects [3,[24][25][26], also confirms the stabilizing effect of refuge in predator-prey systems [35,36].However, the spatiotemporal dynamics of discrete insect predator-prey systems with the refuge effect are still unclear.
In this study, we chose two sets of moth-spider parameters, and then discretized the predator-prey model with a modified Leslie-Gower and Holling type II functional response incorporating the refuge effect [37] to reflect the spatiotemporal dynamics of insect predatorprey systems with relatively lower and higher growth rates.This study is organized as follows: (1) discretize the continuous systems based on coupled image lattice models; (2) apply the center manifold theorem and bifurcation theorems to get the generation conditions of the flip bifurcation and Neimark-Sacker bifurcation; (3) analyze the refuge effect on the spatiotemporal development behavior in moth and spider populations under various types of simulation results, such as bifurcation results, space-time development and spatial amplitude.Our study is intended to reveal the spatiotemporal dynamics of the insect predator-prey system with the refuge effect, which should provide a theoretical basis for future biological control of pests.

Spatiotemporally Discrete Model
In this study, we investigate the modified Leslie-Gower and Holling type II functional response predator-prey model, which incorporates refuge effect with protection for a certain percentage of prey [37]: In which H and P describe prey and predator population densities at time t; r represents the growth rate of prey; a represents the intra-specific competition coefficient of prey; mH represents the role of refuge in the protection of prey, with m representing the proportion of prey protected, mϵ[0, 1] ; H/(b + H) represents the Holling type II functional response; d represents the maximum growth rate of predator; cP/(b + H) represents the modified Leslie-Gower-type numerical response function; c represents the maximum value of the per capita reduction in prey due to the predator; b measures the extent to which the environment provides protection (prey and predator); and r, a, b, c, d are taken as positive constants (see the first table in the Section 2.4).
In reality, refuge comes in all sizes, with smaller refuge only protecting a certain percentage of prey.On this basis, Equation (1) considers the following cases: when m = 0, Model (1) represents there being no the refuge here; when m = 1, Model (1) represents that the refuge here protects all prey; and when 0 < m < 1, Model (1) represents that the refuge protects a certain percentage of prey.
In order to observe richer dynamics in the insect-predator system, we add the space term to Model (1), i.e., In which D 1 and D 2 describe the diffusion coefficient of the predation system and ∇ 2 describes the Laplace operator in two dimensions.
Considering that refuge is spatially distributed in patches, we discretize the model (2).Firstly, the space is divided into N × N grid units, in which t is the time interval and h is the space interval.Each grid represents a spatial patch, which does energy exchange with surrounding patches.By applying the spatiotemporally discrete lattice of coupled images, we derive the spatiotemporally discrete predator-prey system, i.e., where H(i, j, t) and P(i, j, t)(i, j ∈ {1, 2, 3, • • • , n}, t ∈ Z + ) represent the prey and preda- tor densities at time t at the (i, j) lattice point; H ′ (i, j, t) and P ′ (i, j, t) represent the post-diffusion prey and predator densities; τ and δ describe the spatiotemporally discrete time and space intervals; ∇ 2 d is the discrete Laplace operators; and f 1 and f 2 are functions that are determined by intra-and inter-species interactions between prey and predator, which are as follows: The specific description of a spatiotemporally discrete model incorporating the refuge effect is shown by (2)~ (6).In this discrete model, some patches can be randomly selected among N × N patches to indicate refuge of prey, in which we apply Equation (1).This model can reflect some real phenomena and have practical relevance, which lays the groundwork for the parsing analysis for the future steps.

Theorem 1. (a)
The fixed point E 0 and E 1 are unstable under any parameter condition; (b) the fixed point E 2 is not stable when r > (1 − m)d/c; (c) the fixed point E 3 is stable when q < 1 and −(1 Proof.The Jacobian matrix of Map ( 7) can be analyzed for the stability of fixed points [38,39] and its expression at any point is The eigenvalues of the Jacobian matrices J(E 0 ) and J(E 1 ) of the fixed point E 0 and E 1 are λ 1 = 1 + τr and λ 2 = 1 + τr.Since λ 2 > 1, the fixed points E 0 and E 1 are unstable according to Lemma 1.
In addition, the eigenvalues of the Jacobian matrix J(E 2 ) of the fixed point E 2 is Similarly, the eigenvalues of the Jacobian matrix J(E 3 ) of the fixed point E 3 are When q < 1 and −(1 + q) < p < 1 + q, there is |λ 1 | < 1 and |λ 2 | < 1, by Lemma 1 the fixed point E 3 is stable.□ Figure 1 illustrates the parameter space diagram of the fixed points E 1 , E 2 , and E 3 under two sets of parameters.Moreover, from Theorem 1, when E 2 is stable, the discrete predator-prey system will converge to the state where the predator is extinct and prey exists only.The extinction of the predator in this state suggests the predator-prey interactions disappear, which leads to no pattern formation.Therefore, we merely focus on the pattern formation at the stable point E 3 .
According to the center manifold theorem [41,42] and the bifurcation theorems, we conduct a bifurcation analysis of spatiotemporally discrete predator-prey systems around the stable point E 3 using the refuge effect parameter m as the bifurcation parameter.Then, we determine the generate conditions of flip bifurcation and Neimark-Sacker bifurcation, respectively.
Figure 1 illustrates the parameter space diagram of the fixed points  ,  , and  under two sets of parameters.Moreover, from Theorem 1, when  is stable, the discrete predator-prey system will converge to the state where the predator is extinct and prey exists only.The extinction of the predator in this state suggests the predator-prey interactions disappear, which leads to no pattern formation.Therefore, we merely focus on the pattern formation at the stable point  .According to the center manifold theorem [41,42] and the bifurcation theorems, we conduct a bifurcation analysis of spatiotemporally discrete predator-prey systems around the stable point  using the refuge effect parameter m as the bifurcation parameter.Then, we determine the generate conditions of flip bifurcation and Neimark-Sacker bifurcation, respectively.

Flip Bifurcation Analysis
According to the flip bifurcation theorem [41,43], the flip bifurcation occurs when one of the two eigenvalues of ( ) satisfies  = −1.If give  =  − 1, in other words, ( ) ( ) ( ) in which According to the Cardano formula, the above one-dimensional cubic Equation ( 13) can be equated to

Flip Bifurcation Analysis
According to the flip bifurcation theorem [41,43], the flip bifurcation occurs when one of the two eigenvalues of J(E 3 ) satisfies λ = −1.If give q = p − 1, in other words, in which According to the Cardano formula, the above one-dimensional cubic Equation (13) can be equated to (m Then can be obtained as In which With the above conditions (16) satisfied, which leads to the two eigenvalues of J(E 3 ) become λ 1 = −1 and λ 2 = 1 − p(m * ), and the generation of the flip bifurcation requires the satisfaction Entropy 2024, 26, 196 6 of 21 When Equation ( 16) is satisfied, let By translating the fixed point of map (7) to the origin and performing a Taylor expansion near the fixed point E 3 , we can transform map (7) into the following expression: is a polynomial function of at least fourth order for the dependent variable and , Then, according to Map (20), we get a reversible transformation as follows: In which On the basis of the center manifold theorem [41,42], the center prevalence W C (0, 0, 0) of Map (21) exists at the fixed point (0, 0, 0), which can be approximated as In which e 0 = 0, e 1 = 0, e 2 = a 12 (1+a 11 )(a 14 +2a 25 )+(1+a 11 ) 2 (2a 25 −a 28 )+a .(25) Correspondingly, the dynamics of Map (19) restricted to the center manifold W C (0, 0, 0) are obtained, namely: (26)   where , , In addition, the occurrence of a flip bifurcation for Map (20) requires the following two determinants to be nonzero: The calculation gives that Theorem 2. In the discrete predator-prey system (4), a flip bifurcation at the fixed point (H * , P * ) oc- curs if Equations ( 16) and ( 18), µ 2 ̸ = 0 and µ 5 ̸ = µ 2 1 are satisfied.When η 2 > 0, the period-2 points bifurcating from the (H * , P * ) are stable; when η 2 < 0, the bifurcating period-2 points are unstable.
The detailed calculation process of Neimark-Sacker bifurcation is shown below.The parameter values are r = 3.37, a = 0.8, b = 0.9, d = 0.33, and c = 0.44 after calculating the critical value of the occurrence of the flip bifurcation m * = 0.054124.At this time, the discriminant values are η 1 = −0.096and η 2 = 0.4971.On the basis of Theorem 2, we can observe that the period-2 orbitals of the system stabilize from the bifurcation of (H * , P * ).Moreover, from Figure 2a, it can be seen that for the above given parameters, the system is undergoing the flip bifurcation at m * .
In addition, the occurrence of a flip bifurcation for Map (20) requires the following two determinants to be nonzero: The calculation gives that  =  ,  =  +  .
The detailed calculation process of Neimark-Sacker bifurcation is shown below.The parameter values are  = 3.37,  = 0.8,  = 0.9,  = 0.33, and  = 0.44 after calculating the critical value of the occurrence of the flip bifurcation  * = 0.054124 .At this time, the discriminant values are  = −0.096and  = 0.4971.On the basis of Theorem 2, we can observe that the period-2 orbitals of the system stabilize from the bifurcation of ( * ,  * ).Moreover, from Figure 2a, it can be seen that for the above given parameters, the system is undergoing the flip bifurcation at  * .(Around  , the system in all regions has always been unstable and cannot produce bifurcations; around  , the system from the blue region into the red region gradually converges to the state where the predator is extinct; and around  , the system from the blue region into the yellow region is undergoing the flip bifurcation, and from the blue region into the green region is undergoing Neimark-Sacker bifurcation).(Around E 1 , the system in all regions has always been unstable and cannot produce bifurcations; around E 2 , the system from the blue region into the red region gradually converges to the state where the predator is extinct; and around E 3 , the system from the blue region into the yellow region is undergoing the flip bifurcation, and from the blue region into the green region is undergoing Neimark-Sacker bifurcation).According to the Neimark-Sacker bifurcation theorem [44], Neimark-Sacker bifurcation-generating conditions around the spatiotemporally discrete predator-prey first requires that the eigenvalues in Equation ( 10) are a pair of conjugate complexes: Modulo 1, i.e., p 2 − 4q < 0 and q = 1, that is to say, Equation ( 29) is equivalent to In which (31) According to the Caldano formula, the above one-dimensional cubic Equation (30) can be equated to (m Then, the following can be obtained: In which , When the parameter conditions satisfy ( 28) and (33), Map (7) can translate the fixed point (H * , P * ) to the origin through the transformation w = H − H * and z = P − P * .Then, using the Taylor expansion, we can obtain  The two eigenvalues of Map (35) at the fixed point (0, 0) are also conjugate complexes of Mode 1, denoted by where i ' = √ −1, p(m) and q(m) are described in Equations ( 9) and (10).The value in Equation ( 33) is obtained as |λ| = q 2 (m 0 ) = 1 by substituting it into Equation (36).In addition, the condition to generate the Neimark-Sacker bifurcation is also required to be satisfied.
From Condition (38), it follows that Next, in order to study the canonical type of Map (35), we apply the following reversible transformation to Map ( 35) Then, Map (35) transforms to In which In order for Map (41) to undergo a Neimark-Sacker bifurcation near (0, 0), the follow- ing determinant needs to be satisfied that it is not zero: In which Then, bringing Equation ( 46) into Condition (45), we get In which Theorem 3. In the discrete predator-prey system, a Neimark-Sacker bifurcation occurs at the fixed point (H * , P * ) if conditions ( 28), ( 33), ( 37), (39), and (47) are satisfied.When a 0 < 0 and d < 0 , the invariant ring of attraction bifurcates at m < m 0 ; when a 0 > 0 and d > 0, the invariant ring of repulsion bifurcates at m > m 0 .
The detailed calculation process of flip bifurcation is shown below.The parameter values are r = 1.19, a = 0.6, b = 0.16, d = 0.37, and c = 0.42 after calculating the critical value of the occurrence of the Neimark-Sacker bifurcation m 0 = 0.104355.At this time, the discriminant values are d = −0.9219, a 0 = −0.0498,On the basis of Theorem 3, the system undergoes Neimark-Sacker bifurcation from the attractive invariant ring when m < m 0 .Moreover, from Figure 2b, it can be seen that for the above given parameters, the system first experiences Neimark-Sacker bifurcation then enters a steady state at m 0 .

Numerical Simulation Methods
In this study, we perform various numerical simulations (including bifurcation diagrams, maximum Lyapunov exponential diagrams, phase-plane diagrams, spatial amplitude diagrams, and space-time diagrams) to present the results of parsing analysis intuitively and explore the spatiotemporal dynamics of insect predator-prey systems deeply, as the following: Bifurcation diagrams are graphs that show the various bifurcation behaviors of the dynamical system, usually with one parameter as the horizontal coordinate and one state variable as the vertical coordinate.In this study, we use moth and spider population densities as vertical coordinates and the refuge effect parameter as horizontal coordinates.Then, we can graphically get the bifurcation information of population dynamics in response to refuge effects.
The maximum Lyapunov exponent diagram is an important indicator used to determine whether a system is chaotic or not, which quantitatively portrays the average degree of convergence of the distance between two neighboring trajectories in the phase-plane space.For discrete systems H t+1 = f (H t ), the Lyapunov exponent is satisfied δ H t | ≈ e λt |δ H 0 | where δH 0 represents distance between two initial states in the phase-plane space, λ represents the Lyapunov exponent, and δH t represents that the discrete system evolves from two different initial states to two new states at a distance in the phase-plane space after time t.The Lyapunov exponent is given by λ(H 0 ) = lim t→∞ dH ).When the Lyapunov exponent is greater than 0, the system is in a chaotic state at that time.In addition, the higher the Lyapunov exponent, the more obvious the chaotic characteristics and the stronger the degree of chaos.
Furthermore, in order to clearly present the spatial information of the bifurcation map under different degrees of refuge effect, we draw the phase-plane diagram, which reflects the projection of the numerical solution of the dynamical system onto the phase space or its subspace.This allows us to observe the dynamical states of the system, such as fixed points, invariant rings, chaotic attractors, and so on.However, different chaotic phase-plane diagrams are difficult to recognize, so we apply the following methods to differentiate them.
Time series diagrams demonstrate the change of a state variable over time, where the system variable is the average biomass of prey, which shows the stability of the system.Space-amplitude and space-time diagrams are used to characterize the phase-plane diagram and distinguish chaotic phenomena.The space-amplitude diagram is a superposition of the state quantities of each lattice over a period, in which the horizontal coordinate is the spatial location and the vertical coordinate is the population density.The space-time diagram is a plot of the state quantities of each lattice over time, in which the horizontal coordinate is time and its vertical coordinate is spatial location.In this study, the space-amplitude diagram is drawn by superimposing the first 50 spatial lattices after 60,000 iterations and selecting the state quantities from the last 10,000 iterations [44,45].The space-time diagram is drawn by superimposing the amplitudes of 100 spatial lattices and plotting the state quantities every 20 iterations after 50,000 iterations [44,45].From these, we can find the spatiotemporal evolution patterns on chaotic routes of discrete insect predator-prey systems.
All simulations are performed in MATLAB R2016a.By numerical simulation, we observe and discuss the population development dynamics of insect communities.This can then be used to improve pest control and better understand the behavior of insect predatorprey systems.The plotting of these diagrams is based on the theoretical calculations in Section 2.2 and parameter selections in Section 2.4.

Parameter Selection
Moths and spiders form a typical insect interaction, occupying an important position.Moths, due to their destructive power and high fecundity, are well known as destructive pests [46].Spiders, as predators of insects and other invertebrates, are an important component of natural enemies.So, we use the moth-spider parameter to reflect some of the phenomena of insect communities.In this case, the spiders hunt moths by moving around.In addition, in order to more closely simulate the degree of protection for moths in the refuge, we use all space patches as a refuge in one scenario, and randomly select 50% space patches as the refuge in another scenario.The parameter values used in this study were set referring to the prior literature (Table 1).
Subsequently, according to the literature [47][48][49], we selected two sets of parameters for moth populations with relatively higher and lower growth rates, respectively.We then performed bifurcation analyses based on the two sets of parameters, respectively.

Flip Bifurcations Produced by Moths with Relatively Higher Growth Rates
The moth parameter set with higher relative growth rates (see Table 1) produces a flip bifurcation.the flip bifurcation diagrams show that the moth (Figure 3a) and spider (Figure 3b) populations gradually transition from a uniformly fixed state to a chaotic oscillatory state as the parameter m (mϵ(0, 0.45] ) increases.When m < m * , the fixed point (H * , P * ) is stable.However, when m > m * , the phenomenon of period-doubling cascades occurs.From the bifurcation diagrams, it can be observed that the increase in the refuge effect led to a greater population density of moths, while the population density of spiders decreased in a short period.The maximum Lyapunov exponential diagram (Figure 3c) shows that as m increases beyond the critical value of 0.317, the perioddoubling cascades causes the system to exhibit chaotic behavior for the first time.From the maximum Lyapunov exponential diagram, it can be observed that moth and spider population dynamics become chaotic and unpredictable as the refuge effect increases.

Flip Bifurcations Produced by Moths with Relatively Higher Growth Rates
The moth parameter set with higher relative growth rates (see Table 1) produces a flip bifurcation.the flip bifurcation diagrams show that the moth (Figure 3a) and spider (Figure 3b) populations gradually transition from a uniformly fixed state to a chaotic oscillatory state as the parameter  (ϵ(0, 0.45]) increases.When  <  * , the fixed point ( * ,  * ) is stable.However, when  >  * , the phenomenon of period-doubling cascades occurs.From the bifurcation diagrams, it can be observed that the increase in the refuge effect led to a greater population density of moths, while the population density of spiders decreased in a short period.The maximum Lyapunov exponential diagram (Figure 3c) shows that as m increases beyond the critical value of 0.317, the period-doubling cascades causes the system to exhibit chaotic behavior for the first time.From the maximum Lyapunov exponential diagram, it can be observed that moth and spider population dynamics become chaotic and unpredictable as the refuge effect increases.The spatiotemporal behavioral development process of the fixed point, the periodic attractor, and the chaotic attractor generated along the chaotic routes under the flip bifurcation is shown in Figures 4 and 5. Ideally, we selected 100% of the space patches as refuge (Figure 4).The time series diagrams (Figure 4b) show that as  increases, the system shifts from a steady state to the oscillatory and unstable state, eventually becoming irregularly chaotic.The spatial amplitude (Figure 4c) shows that the overall moth The spatiotemporal behavioral development process of the fixed point, the periodic attractor, and the chaotic attractor generated along the chaotic routes under the flip bifurcation is shown in Figures 4 and 5. Ideally, we selected 100% of the space patches as refuge (Figure 4).The time series diagrams (Figure 4b) show that as m increases, the system shifts from a steady state to the oscillatory and unstable state, eventually becoming irregularly chaotic.The spatial amplitude (Figure 4c) shows that the overall moth population in the immobile point state remains constant over time.The system is in a state of cyclic motion at Cycle-2.The orbit of Cycle-4 shows a state of twisting and alternating kinks, in which the positions depend on the initial values.When the system firstly enters the chaotic region (m = m * ), the spatial amplitude diagram starts to show chaotic and synchronized chaotic regions (m = 0.3755), which eventually become difficult to identify.
The space-time diagrams (Figure 4d) in the state of fixed point and periodic attractors show a stable state that is spatially homogeneous, no discontinuities occur, and the pattern remains constant over time.Subsequently, the stable boundary of this system gradually appears fragmentary and defects, and some lattice points even exhibit chaotic behavior.When m increases to 0.3784, the degree of chaos is further strengthened and the spatiotemporal diagrams show multiple stabilized modes transformed by explosive change.It is then difficult to recognize any particular mode in time and space, and the space lattice is almost in a chaotic and turbulent state.The above results indicate that the spatial amplitude and space-time diagrams gradually show chaotic features of disorder and irregularity as the refuge effect increases.In addition, the system gradually develops from a stabilized to fully developed turbulence pattern during this period.It undergoes the frozen random pattern, the defect chaotic diffusion pattern, and the competition intermittency pattern.
More truthfully, we randomly select 50% of the space patches as the refuge (Figure 5).As can be seen in Figure 5a, as m increases, the system gradually becomes unstable and shifts to a chaotic state.Moreover, in comparison to Figure 4c, when the distribution of the refuge is discontinuous, the spatial amplitude changes from a steady state to a state of cyclic motion at the fixed point; the original cyclic motion of the system at Cycle-2 is separated by a number of equilibrium line segments; and the same is true for the Cycle-4 orbital and chaotic regions.The space-time diagrams in the fixed point and periodic attractor states show a wider variety of stable and defect-free frozen bands than in Figure 4d.Under the remaining parameters, it can be seen that new patterns correspond to mixed patterns with frozen bands.
The numerical simulation of flip bifurcation reveals that, as the refuge effect increases to m * , the spatiotemporal dynamics of the moths and spiders shift from a steady state to a chaotic state.Starting from m = 0.317, the structure of prey populations appears chaotic, which increases the persistence of population survival.These are consistent with previous results "an increase in refuges can increase the density of prey populations, as well as trigger large-scale population explosions [36,51,52]".In specific conditions, a higher refuge effect contributes to the reproduction of moths.But it may reduce spider population density in a short period, which is consistent with empirical studies [53,54].Furthermore, through the time series diagrams, the population density of moths gradually tends to stabilize with the development of time when the refuge effect is small.As the refuge effect increases, the population density of moths becomes more irregular and chaotic over time.This suggests that increasing the refuge effect under certain conditions decreases the stability of the system, which is verified by the numerical simulations of the previous theoretical results [55].Advances on previous studies are made by analyzing the spatial amplitude and space-time diagrams of the system [45,55,56].We further refine four spatiotemporal development phases (from the frozen random pattern to the defect chaotic diffusion pattern, then the competition intermittency pattern, and finally to the fully developed turbulence pattern) in the chaotic regions and give a means to study spatiotemporal behavior in chaotic population states.Moreover, as the number of refuges decreased, we also found new patterns of the frozen random pattern mixed with the defect chaotic diffusion pattern, competition intermittency pattern, and fully developed turbulence pattern, respectively.
Moth populations become stable at regular intervals as the number of refuges decreases under the same refuge effect.These findings offer theoretical guidance (e.g., removing their refuges or setting traps on their escape routes in time) for controlling the reproduction of pests with abundant food sources and high growth rates.The numerical simulation of flip bifurcation reveals that, as the refuge effect increases to  * , the spatiotemporal dynamics of the moths and spiders shift from a steady state to a chaotic state.Starting from  = 0.317, the structure of prey populations appears chaotic,

Neimark-Sacker Bifurcation Produced by Moths with Relatively Lower Growth Rates
The moth parameter set with relatively lower growth rates (see Table 1) produces a Neimark-Sacker bifurcation showing that the moth (Figure 6a) and spider (Figure 6b) populations change in the interval of mϵ(0, 0.15] .When m > m 0 , the fixed point (H * , P * ) is stable.However, the Neimark-Sacker bifurcation occurs when m < m 0 , leads to the destabilization of the moth-spider system.We find that as the refuge effect increases, moths with relatively lower growth rates remain stable at relatively low densities.Moreover, the maximum Lyapunov exponential diagram (Figure 6c) demonstrates that the system stays in an unchaotic state after the Neimark-Sacker destabilization.
and fully developed turbulence pattern, respectively.Moth populations become stable at regular intervals as the number of refuges decreases under the same refuge effect.These findings offer theoretical guidance (e.g., removing their refuges or setting traps on their escape routes in time) for controlling the reproduction of pests with abundant food sources and high growth rates.

Neimark-Sacker Bifurcation Produced by Moths with Relatively Lower Growth Rates
The moth parameter set with relatively lower growth rates (see Table 1) produces a Neimark-Sacker bifurcation showing that the moth (Figure 6a) and spider (Figure 6b) populations change in the interval of ϵ(0, 0.15].When  >  , the fixed point ( * ,  * ) is stable.However, the Neimark-Sacker bifurcation occurs when  <  , leads to the destabilization of the moth-spider system.We find that as the refuge effect increases, moths with relatively lower growth rates remain stable at relatively low densities.Moreover, the maximum Lyapunov exponential diagram (Figure 6c) demonstrates that the system stays in an unchaotic state after the Neimark-Sacker destabilization.After the destabilization of the Neimark-Sacker bifurcation, the phase-plane diagrams show that the moth-spider system only exhibits spatiotemporal behavior characterized by the fixed point and the invariant ring (Figure 7a).Ideally, we selected 100% of the space patches as the refuge (Figure 7).As the value of  decreases, the system goes from a stable state to a regularly unstable state (Figure 7b) and undergoes only a frozen random pattern throughout.The corresponding space-amplitude variations display a regular clustered structure.The space-time diagrams show that stable band and After the destabilization of the Neimark-Sacker bifurcation, the phase-plane diagrams show that the moth-spider system only exhibits spatiotemporal behavior characterized by the fixed point and the invariant ring (Figure 7a).Ideally, we selected 100% of the space patches as the refuge (Figure 7).As the value of m decreases, the system goes from a stable state to a regularly unstable state (Figure 7b) and undergoes only a frozen random pattern throughout.The corresponding space-amplitude variations display a regular clustered structure.The space-time diagrams show that stable band and bar structures with clear and frozen boundaries do not vary with time.Therefore, the spatiotemporally discrete moth-spider system maintains relative stability both in time and space when the growth rates of moths are relatively lower.
More truthfully, we randomly select 50% of the space patches as the refuge (Figure 8).The time series diagrams in Figure 8a show that as m decreases, the system is regularly but not chaotically unstable.Moreover, compared to Figure 7c, the spatial amplitudes are separated by a number of equilibrium line segments, showing a block of knotty organization.The space-time diagrams show intermittent banding and bar structures with clear and frozen boundaries do not vary with time.
The numerical simulation of Neimark-Sacker bifurcation reveals that the population of moths and spiders progressively transforms from oscillatory behavior to stability with the refuge effect increasing, suggesting that the refuge effect exerts a stabilizing influence on both moth and spider populations, which is in agreement with most previous results [51,55].The population pattern of the moth system at low growth rates remains in a regularly cycling spatiotemporal state.In addition, as the number of refuge decreases, moth populations become stable at regular intervals under the same refuge effect.Our results complement previous conclusions "increasing refuge effects can trigger population outbreaks" [36,51,52].These provide theoretical guidance for the management of insect systems.For beneficial insects with a lower growth rate and a preference for specific food sources, appropriately managing the refuge effect is beneficial to maintaining their population stability.However, for pests with low growth rates, we should remove the refuge at times to destabilize their populations.
bar structures with clear and frozen boundaries do not vary with time.Therefore, the spatiotemporally discrete moth-spider system maintains relative stability both in time and space when the growth rates of moths are relatively lower.More truthfully, we randomly select 50% of the space patches as the refuge (Figure 8).The time series diagrams in Figure 8a show that as  decreases, the system is regularly but not chaotically unstable.Moreover, compared to Figure 7c, the spatial amplitudes are separated by a number of equilibrium line segments, showing a block of knotty organization.The space-time diagrams show intermittent banding and bar structures with clear and frozen boundaries do not vary with time.The numerical simulation of Neimark-Sacker bifurcation reveals that the population of moths and spiders progressively transforms from oscillatory behavior to stability with the refuge effect increasing, suggesting that the refuge effect exerts a stabilizing influence on both moth and spider populations, which is in agreement with most previous results [51,55].The population pattern of the moth system at low growth rates remains in a

Conclusions
Insect predator-prey systems are extremely easy to produce rich dynamic behaviors with, of which the moths-spider system is a typical example.In this study, we clearly unveil the complex spatiotemporal behavior and development law of the insect populations incorporating a refuge effect.Via numerical simulation of two sets of parameters in the moth-spider interaction, we revealed the dynamical behaviors of moths with different growth rates under different levels of the refuge effect and further reflecting the complex mechanisms of insect predator-prey systems, which are shown as follows: For the moth parameter group with relatively high growth rates, the system produces a flip bifurcation and opens the routes towards chaos.The flip bifurcation appears in the period-2, 4, 8, and 16 bifurcation cascades, a periodic window of period-6 orbits, and chaotic sets.Along these chaotic paths, the system exhibits four distinct spatiotemporal behavior phases: the frozen random pattern, the defective chaotic diffusion pattern, the competition intermittency pattern, and fully developed turbulence pattern.As the refuge effect increases, the initially stable frozen spatiotemporal boundaries break down, then develop into an overall ordered but locally chaotic competitive state, and eventually to chaotic turbulence.And under certain conditions increasing the refuge effect decreases the stability of the system.Moreover, as the number of refuges decreases, we also found new patterns of frozen random pattern mixed with defect chaotic diffusion pattern, competition intermittency pattern, and fully developed turbulence pattern, respectively.This means that the refuge effect can increase the persistence of insect population survival and induce disorganization in spatial and temporal patterns.
For the moth parameter group with relatively low growth rates, the system produces an inverse Neimark-Sacker bifurcation, while it does not produce chaos.This system only appears in the form of fixed points and invariant rings, where the spatiotemporal dynamics always show the frozen random pattern.As the refuge effect increased, moths are able to remain stable at low densities.In addition, as the number of refuge decreases, stabilized space-time structures become more fragmented.This implies that properly reducing the refuge effect can destabilize insect populations, and insect systems under such parameters are always spatially and temporally stable.
In nature, insects often employ various strategies to evade their natural enemy, such as seeking refuges in rocks and bushes, or escaping from predator in time through effective information transmission among insects.In this study, the model explored has practical reference values in insect communities under the refuge effect, and numerical simulations (moths with different growth rates) can yield theoretical support for adjusting refuge effect to destabilize pest populations.Our study suggests that regulating the refuge effect can effectively control pest density, which provides valuable insights for managing invertebrate population dynamics or insect pollination.

Figure 2 .
Figure 2. (a) Parameter space diagram (, ) of  ,  , and  under  = 0.8,  = 0.9,  = 0.44,  = 0.33, (b) parameter space diagram (, ) of  ,  , and  under  = 0.6,  = 0.16,  = 0.42,  = 0.37.(Around , the system in all regions has always been unstable and cannot produce bifurcations; around  , the system from the blue region into the red region gradually converges to the state where the predator is extinct; and around  , the system from the blue region into the yellow region is undergoing the flip bifurcation, and from the blue region into the green region is undergoing Neimark-Sacker bifurcation).

Figure 2 .
Figure 2. (a) Parameter space diagram (m, r) of E 1 , E 2 , and E 3 under a = 0.8, b = 0.9, c = 0.44, d = 0.33, (b) parameter space diagram (m, r) of E 1 , E 2 , and E 3 under a = 0.6, b = 0.16, c = 0.42, d = 0.37(Around E 1 , the system in all regions has always been unstable and cannot produce bifurcations; around E 2 , the system from the blue region into the red region gradually converges to the state where the predator is extinct; and around E 3 , the system from the blue region into the yellow region is undergoing the flip bifurcation, and from the blue region into the green region is undergoing Neimark-Sacker bifurcation).
which the coefficients are all given in Equation (20) with replacements.

Figure 7 .
Figure 7. (a) Phase-plane diagrams, (b) time series diagrams, (c) space-amplitude diagrams, and (d) space-time diagrams.The parametric conditions are the same as those in Figure 3, except  = 0.03, 0.1, and 0.11, respectively, from the first to the last row.(selecting of 100% of space patches as the refuge).

Figure 7 .
Figure 7. (a) Phase-plane diagrams, (b) time series diagrams, (c) space-amplitude diagrams, and (d) space-time diagrams.The parametric conditions are the same as those in Figure 3, except m = 0.03, 0.1, and 0.11, respectively, from the first to the last row (selecting of 100% of space patches as the refuge).

Figure 8 .
Figure 8.(a) Time series diagrams, (b) space-amplitude diagrams, and (c) space-time diagrams.The parametric conditions are the same with that in Figure 3, except  = 0.03, 0.1, and 0.11, respectively, from the first to the last row.(selecting of 50% of space patches as refuge).

Figure 8 .
Figure 8.(a) Time series diagrams, (b) space-amplitude diagrams, and (c) space-time diagrams.The parametric conditions are the same with that in Figure 3, except m = 0.03, 0.1, and 0.11, respectively, from the first to the last row (selecting of 50% of space patches as refuge).