Time delays shape the eco-evolutionary dynamics of cooperation

We study the intricate interplay between ecological and evolutionary processes through the lens of the prisoner’s dilemma game. But while previous studies on cooperation amongst selfish individuals often assume instantaneous interactions, we take into consideration delays to investigate how these might affect the causes underlying prosocial behavior. Through analytical calculations and numerical simulations, we demonstrate that delays can lead to oscillations, and by incorporating also the ecological variable of altruistic free space and the evolutionary strategy of punishment, we explore how these factors impact population and community dynamics. Depending on the parameter values and the initial fraction of each strategy, the studied eco-evolutionary model can mimic a cyclic dominance system and even exhibit chaotic behavior, thereby highlighting the importance of complex dynamics for the effective management and conservation of ecological communities. Our research thus contributes to the broader understanding of group decision-making and the emergence of moral behavior in multidimensional social systems.

www.nature.com/scientificreports/ in which the entries represent the payoff accumulated by the player in the left.We consider an additional strategy, the 'punisher' who acts like a cooperator and receives the same payoff R = 1 if the other player decides to cooperate.However, if the opponent player defects, the punisher will use their own resources to punish them and earn a payoff value of S − δ = −δ .In return, the defector makes T − δ = β − δ with δ > 0 .Note that a defector earns more if they interact with a cooperator.This inclusion of punishers will extend our payoff matrix (1) to the following matrix We further incorporate the ecological contribution of free space in the payoff matrix.Free space allows others to replicate and provide others ample opportunity to survive.Interestingly, free space never anticipates anything in return, and this selfless contribution of free space motivates us to update our payoff matrix (2) in the following way, From the above upgraded payoff matrix (3) and in evolutionary sense, ecological free space (F) seems out to be one selfless strategy to be made by any player population and on interaction with this, each population kind gets positive attributes as their benefits for replication, whereas, free space stands out to exist independently for every player kinds to shower it's altruistic ease, without being chosen by any trait as their action in evolution.In order to showcase our dynamical system, we use the positive attributes as benefits to the competing players from this spatial free space but in game theoretic sense, free space does not act as an evolutionary move.Here, all these parameters σ 1 , σ 2 , and σ 3 are strictly positive quantities, as free space contributes altruistically to all the rational individuals.Free space earns only zero in return as it never expects any benefits for its generous acts.Now we calculate each strategy's fitness and determine the evolution of populations by assuming that an individual's reproduction rate depends solely on their average fitness.Let us assume x, y, z, and w be the respective fraction of cooperators, punishers, defectors, and free space, respectively.Hence, x, y, z, w ∈ [0, 1] and x + y + z + w = 1 .Interestingly, the process of learning takes time and effort.It takes a considerable amount of time for people to adapt strategies based on the information they learn.The players need to gather the information at each round of the game and then assess the effectiveness of the methods.Based on their understanding, they spread that information which helps others to recognize which strategy is the most successful in society.Thus individuals select a strategy at time t based on the fitness before τ ≥ 0 time instance.We introduce the variables x τ = x(t − τ ) , y τ = y(t − τ ) , z τ = z(t − τ ) and w τ = w(t − τ ) .Since we are interested in inferring the obtained results using fundamental principles of biological systems, we maintain the constraint x τ + y τ + z τ + w τ = 1 along with x τ , y τ , z τ , w τ ∈ [0, 1] throughout the article.This constraint helps us to eliminate one independent variable and construct a simple eco-evolutionary model.The respective average fitness of cooperators, punishers, defectors and free space is given by We assume all players die with a uniform rate ξ > 0 .Thus, our proposed delayed system looks like Thus substituting Eq. (4) in Eq. (5), we obtain Note that τ = 0 gives the non-delayed system.In the subsequent section, we comprehensively discuss the difference between the outcomes in delayed and non-delayed systems. (2) (3) (4)

Numerical results
Comparison between delayed and non-delayed model.To understand the impact of time delay in population dynamics, specifically in the context of a prisoner's dilemma game, we present the comparative temporal behavior of variables in the presence and absence of delay in Fig. 1.To start with, we consider the values of the parameters as ξ = 0.50 , β = 2.15 , δ = 1.4 , σ 1 = 0.52 , σ 2 = 0.72 and σ 3 = 0.41 in the subfigures (a,b).Subfig- ures of (a) demonstrate small amplitude oscillations of all variables without delay.Clearly, the chosen parameter values help the defectors to dominate others as we identify the inequality z > x > y .It should be noted that our observation may alter for a different set of initial conditions, as the system (6) is multistable.Later, we will examine the role of initial conditions in our model in detail.We fix the initial conditions at (x 0 , y 0 , z 0 ) = (0.3, 0.3, 0.3) .We further set the initial conditions for the delay variables as x τ (0) = 0.25, y τ (0) = 0.25 , and z τ (0) = 0.25 for the subfigures in the presence of delay.The chosen parameter values and initial conditions provide an opportunity to maintain biodiversity by allowing the coexistence of all strategies in subfigures (a,b).The range of oscillations for x ∈ [0.0775, 0.0779] , y ∈ [0.004732, 0.00477] , and z ∈ [0.1145, 0.1148] is petite in subfigures (a), where the sum of the population lies within the range [0.197, 0.1975].Since the punishment parameter's value (δ = 1.4) is high enough, punishers can not afford to survive in the long run, despite the free space-induced benefits towards the punishers being higher as per our chosen parameter values (σ 2 > σ 1 > σ 3 ) .Since the temptation to defect is high for the defectors as β = 2.15 , defectors are able to overcome the hurdle in the long run in our model.We further analytically calculate the interior equilibrium point (0.0777, 0.0047, 0.1147)) for this set of parameter values, which is found to be an unstable focus node as the eigenvalues of the Jacobian at this point are 1 = 0.1236 , and 2,3 = −0.004± 0.0785i where i = √ −1 .To investigate the delay effect, we allow a small amount of delay τ = 0.24 in the state variables in subfigures (b) by keeping fixed all the parameters' values and initial condi- tions of subfigures (a).This inclusion of delay will not alter the inequality z > x > y observed in subfigures (a); Figure 1.A comparison of time series between delay-free and delayed systems: The time series in this figure compares the dynamics of the system (6) with and without delay.(a) and (c) Correspond to the dynamics of the non-delayed system for τ = 0 , that displays either a small-amplitude oscillation in its three variables' temporal evolution (a) or convergence to the interior equilibrium (c) owing to the choice of different parameter values.In contrast, the inclusion of delay generates larger oscillations in both cases, corresponding time series are displayed in (b) and (d) for chosen values of the delay parameter τ = 0.24 and τ = 1.5 , respectively.Temporal evolution of population fraction of the cooperators, punishers and defectors are depicted in the first, second and third columns, respectively, while the evolution of the total population ( x + y + z ) is shown in the last column.Non-zero delay may trigger the total population fraction to cross the upper bound of unity, leading to an overcrowded solution as shown in the last panel.The coexistence of all three populations either in oscillation or in equilibrium helps maintaining the biodiversity.The initial conditions for the non-delayed variables are set at (0.3, 0.3, 0.3), while those for the delayed state variables are set at (0.25, 0.25, 0.25) in all subfigures.The parameter values are fixed at ξ = 0.50 , β = 2.15 , δ = 1.4 , σ 1 = 0.52 , σ 2 = 0.72 and σ 3 = 0.41 for panels (a-b), and ξ = 0.80 , β = 1.25 , δ = 0.30 , σ 1 = 1.00 , σ 2 = 1.50 and σ 3 = 0.60 for panels (c,d).
A striking difference is also observed for the subfigures (c,d) with the fixed parameter values ξ = 0.80 , β = 1.25 , δ = 0.30 , σ 1 = 1.00 , σ 2 = 1.50 , and σ 3 = 0.60 and fixed initial conditions x 0 = 0.3 , y 0 = 0.3 and, z 0 = 0.3 .The subfigure is drawn additionally with the delay parameter τ = 1.5 with fixed initial condition x τ (0) = 0.25, y τ (0) = 0.25 , and z τ (0) = 0.25 .Subfigure (c) reveals the coexistence of all strategies, and we ana- lytically calculate the interior equilibrium (0.274, 0.407, 0.2) for the chosen parameter values.We find this steady state is a locally stable focus node as the eigenvalues of the Jacobian at this steady state are 1 = −0.0208 ,2,3 = −0.3327± 0.2251i .Interestingly, the punishers' population is the dominant in this case as we identify the inequality y > x > z .The lower abundance of defectors is due to the choice of insufficient temptation to defect (β = 1.25) and free space-induced benefits towards defectors σ 3 = 0.60 is lower.Punishers overcome the fierce struggle as punishment parameter value δ = 0.30 is chosen sufficiently low.Punishers use their own resources to penalize the defectors, and since δ is chosen low here, punishers' resources are not overly utilized.Further- more, our selected parameter values suggest free space-induced benefits toward punishers are higher compared to others (σ 2 > σ 1 > σ 3 ) .The addition of a time delay of suitable strength not only destroys the stability of this interior point but also yields an overcrowded solution as x + y + z exceeds 1.This introduction of time delay leads to an instantaneous change in the system's dynamics.
The emergence of these oscillations hints us the spontaneous emergence of cyclic dominance among those strategies.In the following subsection, we discuss this occurrence in more detail.

Cyclic dominance.
We consider the same parameter values ξ = 0.50 , β = 2.5 , δ = 1.39 , σ 1 = 0.52 , σ 2 = 0.72 , σ 3 = 0.41 and τ = 0.019 and integrate our model ( 6) using Huen's method 79 with 10 7 number of iterations.After discarding a sufficiently long transient of length 9.9 × 10 6 iterations, we present their delayed eco-evolutionary dynamics in Fig. 2. We maintain the same initial conditions as in Fig. 1b.Interestingly, the delay parameter helps the system maintain oscillatory dynamics that attest to the emergence of cyclic dominance in our model.Cyclic dominance 80 allows each strategy to dominate others for a specific time window, which is impossible if we attain steady-state dynamics.The periodic dynamics of each variable portray the coexistence of competing strategies.The maxima of each variable in the first row in Fig. 2 unveils z max > x max > y max .Since our chosen value of the temptation parameter, β = 2.50 is higher than the payoffs of a cooperator and punisher, z max can attain such larger values in its temporal dynamics.On the other hand, the punishment parameter, δ , is fixed at 1.39, which helps to control the defection by penalizing them; however, the punisher also loses this amount while punishing the defectors.This makes punishers vulnerable in society for our specific choices of initial conditions and parameter values.Furthermore, the interplay between all parameters is more pronounced in this figure, as despite altruistic free space contributing more towards the punishers and lesser towards the defectors (σ 2 > σ 1 > σ 3 ) , punishers are still unable to take over the defectors.The temptation to defect is so high that it helps to ignore the selfless contribution of free space, and thus defectors can gain a higher density in the long run.Nevertheless, periodic oscillation allows a window of opportunity for the cooperators to invade the defectors, who in turn can invade the punishers and, ultimately, punishers are able to overrun the cooperators.In this way, cyclic dominance emerges spontaneously and captures the beauty of governing eco-evolutionary dynamics.We further confirm that the obtained solution and their sum remain bounded within [0, 1], providing a biologically feasible solution.The two and three-dimensional projections of the attractor, along with the direction of the flow, are plotted in the second row of Fig. 2.
In the next subsection, we will scrutinize the underlying mechanism that drives the system toward an oscillatory behavior.
Exploring the role of delay in the eco-evolutionary dynamics.To further study the effect of delay in our model (6), we keep all parameter values fixed at ξ = 1.2 , β = 1.6 , δ = 0.30 , σ 1 = 1.35 , σ 2 = 1.50 , σ 3 = 1.35 and vary the delay parameter τ within the closed interval [8,12] with a fixed step length 0.005 and plot the bifurcation diagrams in Fig. 3.We find out the cooperators extinct in the long run throughout the interval, as x remains at zero in Fig. 3b.Nevertheless, the dynamics of the punishers and defectors offer a great variety of dynamical behaviors.Our system experiences a series of transitions from a periodic to a chaotic state as τ is varied in Fig. 3.All these figures provide valuable insight into how our system responds to changes in the delay parameter.
We identify punishers overrule the defectors in the steady state regime for the chosen initial conditions and parameter values.Assuming an initial condition of (0.1, 0.2, 0.5) for the state variables and (0.3, 0.3, 0.3) for the delayed variables, we performed 10 7 iterations of the system (6).We discarded the first 9.8 × 10 6 iterations to ensure that our analysis focuses only on the system's long-term behavior.However, this steady state loses its stability and gives rise to a periodic solution beyond a specific value of τ .We will provide a thorough analysis later to detect the point of this Hopf bifurcation analytically.Interestingly, the range of y is vast compared to that of z (c.f.subfigures (c,d) of Fig. 3).This indicates the punishers gain some kind of opportunity through our setup and can dominate the defectors.Notably, the free space-induced benefits towards the punishers are slightly larger than others (σ 2 > σ 1 = σ 3 ) in this figure.Interestingly, we observe this periodic solution loses its stability and gives rise to a new periodic solution with twice the period of the original one as the delay parameter value is increased.As τ is further increased, the system goes through additional period-doubling bifurcations, giving rise to periodic solutions with four times, eight times, and so on, the period of the original periodic attractor.Eventually, the system enters a chaotic regime, where the dynamics are unpredictable and sensitive to small perturbations.To further validate our findings, we use the Lyapunov exponent to measure the sensitivity of our eco-evolutionary model (6) to small perturbations in its initial conditions.The calculation of Lyapunov exponents for delayed www.nature.com/scientificreports/systems is generally more complicated than that of non-delayed systems due to the need to consider the effect of time delay on the system's dynamics.We calculate the Lyapunov exponents of the system for τ ∈ [8,12] by increasing τ with a fixed step length 0.04.Since our system (6) contains only one delay parameter, the number of Lyapunov exponents is equal to the dimensionality of the non-delayed system (i.e., the number of variables in the system).We plot each of these Lyapunov exponents in Fig. 3a.The first and the second largest Lyapunov exponents reveal valuable information about our system.The largest Lyapunov exponent 1 , shown in deep purple, remains negative up to a certain value of τ , indicating stable steady states.For a suitable range of the delay parameter τ , it converges to zero, which reflects a standard signature of having a periodic solution in our system (6).In both of these stable states, nearby trajectories converge towards each other, meaning that small perturbations in the system's initial conditions lead to similar outcomes.Beyond this range of τ , it offers only the positive value describing the rate of exponential divergence of nearby trajectories.This positive Lyapunov exponent 1 identifies the onset of chaos and validates our bifurcation diagram, which reflects the period doubling route to chaos.Interestingly, most of the previous studies on Lyapunov exponents only concentrate on the largest Lyapunov exponent and ignore the possibility of extracting invaluable information from the other Lyapunov exponents.The second-largest Lyapunov exponent 2 , shown in cyan color dashed line in Fig. 3a, discloses fascinating details on the bifurcation point in our study.It initially remains at negative values when the system's behavior remains constant over time.As the system experiences Hopf bifurcation, it attains a value of zero and again provides only negative values for a range of τ .At the Hopf bifurcation point, the system changes from having a stable steady state to a stable periodic orbit, as observed in y and z variables in subfigures (c,d) of Fig. 3.The second-largest Lyapunov exponent 2 again takes the zero value while the system Figure 2. Cyclic dominance among three competing strategies: (a) Temporal dynamics of the oscillatory coexistence of all three strategies exhibiting the emergence of cyclic dominance among the species that allows each speices to dominate others.The following parameter values are chosen for numerical simulation: ξ = 0.5 , β = 2.5 , δ = 1.39 , σ 1 = 0.52 , σ 2 = 0.72 , σ 3 = 0.41 and, τ = 0.019 .Higher value of the temptation parameter β facilitates the defectors although the punishment parameter δ helps controlling their population by penalizing them.On the other hand cooperators get more benefit from the altruistic free space compared to the defectors due to our chosen parameter values, which helps in maintaining the coexistence of all three strategies.Phase portraits of the periodic attractor along with the direction of the flow in (b-d) two dimensions and in (e) three dimension.The initial values for the state variables associated with the player strategies are considered to be, (x 0 , y 0 , z 0 ) = (0.3, 0.3, 0.3) and those for the delayed variables, are considered to be, (x τ (0), y τ (0), z τ (0)) = (0.25, 0.25, 0.25).www.nature.com/scientificreports/undergoes a period-doubling bifurcation.Furthermore, it acquires the value of zero during the appearance of a new stable periodic solution with twice the period of the original periodic attractor.In this way, the secondlargest Lyapunov exponent helps to validate the points where period-doubling bifurcation occurs.It assumes a value of zero when τ is larger and close to 12, as the system experiences a series of period-doubling bifurcations.The third Lyapunov exponent 3 shown in light pink dotted line in Fig. 3a always remains negative throughout the interval of investigation.All these Lyapunov exponents help characterize the system's dynamics and provide a better understanding of its sensitivity to initial conditions and predictability.
In spite of obtaining such captivating dynamics, we can not infer a few of these results from the biological viewpoint.The dynamics must be bounded within [0, 1]; otherwise when the total population exceeds the maximum value 1, we obtain an overcrowded solution.We plot a dash-dotted horizontal line in subfigures (a) and (c) of Fig. 3, below which the dynamics remain bounded within [0, 1].Within the feasible range, we observe steady states as well as periodic oscillations in this figure.Apart from this dynamical behavior, our system (6) may depict quasiperiodic oscillation too, under a suitable choice of parameter values and initial conditions.The temptation parameter always plays a crucial role in determining the fate of defectors.If it is very large, defectors gain massive benefits that are too large to overcome.Under that circumstances, defectors are the unbeatable winners.If the temptation parameter is moderate, one can think about any approach to overcome the natural selfish instinct.To understand the effect of the temptation parameter in our study, we keep fixed all the parameter values and initial conditions as in Fig. 3 and choose three different pairs of values of (β, τ ) .The dynamics offer a plethora of new behavior depending on these choices.We present all these observations in Fig. 4. Clearly, the temptation to defect can not help cooperators survive and evolve.Thus, we are never able to observe the survival of cooperators in Fig. 4 for the chosen parameter values and initial conditions.The cooperators' fraction x always remains at zero.However, the dynamics of y and z again offer a wide variety.The chaotic dynamics observed in The effect of the delay parameter τ on the dynamics of the total population x + y + z of all three species, that indicates the emergence of periodic oscillation via Hopf bifurcation at τ = 8.99 which then leads to period doubling bifurcation as delay increases and eventually to chaos.The three Lyapunov exponents 1,2,3 of the delayed system are also displayed among which 1 (solid purple) is the maximum that validates whether the dynamics is a steady state, periodic or chaotic oscillation contingent to its value being negative, zero or positive, respectively.The second largest exponent 2 (cyan dashed) remains negative throughout the delay interval except at the bifurcation points where it takes zero value and 3 (pink dotted) remains negative throughout the interval.The individual bifurcation diagrams are also shown in (b) for cooperators x, (c) for punishers y, and (d) for defectors z, respectively.The following parameter values are chosen: ξ = 1.20, β = 1.60, δ = 0.30, σ 1 = 1.35, σ 2 = 1.50, σ 3 = 1.35 and the initial conditions are set at (0.1, 0.2, 0.5) for the non-delayed variables and (0.3, 0.3, 0.3) for the delyaed variables.The figures are plotted by taking the τ step length of 0.005 for the bifurcation diagrams and 0.04 for the Lyapunov exponents in the range [8,12].From the individual bifurcation diagrams it is clear that the chosen parameter values correspond to the coexistence of only the punishers and defectors, with punishers being the dominating population.Fig. 3 is an overcrowded solution as x + y + z > 1 ; hence we can not offer any biological interpretations for such an attractor.The upper panel of Fig. 4 shows a chaotic oscillation in y and z for β = 1.93 and τ = 11.64 ; where x + y + z ∈ [0, 1] .We plot the power spectrum corresponding to this temporal dynamics, which exhibits a broad range of frequencies in their power spectra, with no single dominant frequency.To avoid repetition, we chose not to display the Lyapunov exponent plot in this instance, but we have confirmed the validity of our results through its inclusion.Here for this set of parameter values, there is no ultimate winner among punishers and defectors, as they dominate one another in an unordered way, as revealed in subfigure (d).The chaotic oscillation never allows strict dominance over one another.Nevertheless, this scenario may alter for a different choice of (β, τ ) , and one may dominate another in the long run.We choose β = 1.64 and τ = 10.73 for the middle panel of this figure.The change in two parameter values does not alter the fortune of the cooperators, and they remain extinct in the long run, as shown in subfigure (f).However, maximum values of punishers' density can dominate the defectors' density for this set of parameter values (See subfigures (g-i)).We identify the oscillations that portray two-periodic dynamics.We plot the power spectrum of the time series in subfigure (j) and find these dominant frequencies correspond to the frequencies at which the signal repeats itself.All these subfigures confirm that the system (6) spends some time oscillating with one amplitude or frequency, then switching to the other amplitude or frequency, and continue oscillating in this way for these parameter values and initial conditions.The phase portrait in the yz plane in subfigure (i) of Fig. 4 also confirms this behavior.
The dynamical behavior is much more complex in the bottom panel of Fig. 4 where we choose β = 1.9 and τ = 10.5 .Cooperators are still unable to overcome the fierce struggle for resources and end up saturating in zero density, as shown in Fig. 4k.Punishers and defectors keep oscillating at two or more incommensurate frequencies, meaning that their ratio is irrational.The oscillation frequencies of these two populations exhibit a pattern that repeats itself over time, but the repetition is not strictly periodic (See subfigures (l-n)).The power spectrum of www.nature.com/scientificreports/this attractor in of Fig. 4o confirms the appearance of a broad band of frequency components rather than distinct peaks at specific frequencies.And the width of these bands indicates the degree of incommensurability between the frequencies.This non-periodic and non-chaotic temporal behavior is a subject of deeper investigation and can provide insights into the behavior of complex systems.We would like to investigate the underlying mechanisms driving the system's quasiperiodic dynamics and identify the key factors influencing its behavior in the future.Figure 4 highlights the importance of these two parameters β and τ .In the following subsection, we will investigate the interplay of these two parameters in much more detail.
Investigating the complex relationship between the temptation parameter and delay parameter.Our proposed eco-evolutionary delay model explores various scenarios and player population diversities in society by changing the parameters in the model.The time delay parameter, τ , is critical in altering the population dynamics, while the payoff parameters determine the viability of competing strategies in society.We analyze the survivability of different strategies and their dynamical behaviors by simultaneously varying two parameters: the temptation payoff parameter, β , and the time delay parameter, τ , in a two-dimensional param- eter space.We investigate the behavior of our model ( 6) by varying β ∈ (1, 2] and τ ∈ [8,12] with fixed step length 0.01, while keeping other parameters fixed at ξ = 1.20 , δ = 0.30 , σ 1 = 1.35 , σ 2 = 1.50 , and σ 3 = 1.35 .The initial fractions for the non-delayed state variables are (0.1, 0.2, 0.5), while the initial fractions of the delayed variables are (0.3, 0.3, 0.3).This two-dimensional parameter space in Fig. 5 provides plenty of information about our proposed delayed system (6).In the classical prisoner's dilemma game, defectors have a primary advantage over cooperators.We consider the chosen advantages from the free space as attributes to be equal for both the cooperator and defector populations ( σ 1 = σ 3 = 1.35 ).Moreover, the initial fraction of the population of coop- erators is much lesser than others.Consequently, we end up with a society free from cooperators for any choice of temptation and delay parameters.
Since punishers receive significantly higher benefits ( σ 2 = 1.50 ) from the free space compared to others, only the punishers are observed to survive periodically (indicated by mauve region (A) and by the time series presentation in sub-figure (a) in Fig. 5) while cooperators and defectors can not compete.In this regard, the fraction of the punisher population exceeds the threshold 1.When we vary the time delay by keeping the value of β fixed, we observe an increase in the amplitude of oscillation but no change in the behavior of the competing strategies.However, an increment of β always helps the defectors back in the contest.We locate a small area highlighted by lime yellow (B) in Fig. 5 where defectors can coexist with the punishers.Beyond β = 1.4 , the punishers and defectors begin to dominate one another depending on the time window and yield a periodic attractor here.Sub-figure (b) of Fig. 5 represents the periodic overcrowded oscillations among punishers and defectors.This periodic oscillation occurs due to the Hopf bifurcation, and the value of τ , where the Hopf bifurcation occurs, gradually decreases as the value of β increases.Notably, Ref. 13 also found this critical value of β = 1.40 , beyond which defectors are very hard to defeat.As we further increase the temptation to defect, we observe an attractor with period-2 emerges at β = 1.46 (See the set of sub-figures (j) along with the deep grey area (J) of Fig. 5).We further add the temporal behavior of different regions of the parameter space in Fig. 5.We observe the z max increases as β increases.In other words, the progressive growth of β facilitates the growth of the defector popula- tions.However, the solution remains overcrowded as y + z exceeds the upper bound of unity.
Noticeably, we are able to detect a fair portion, shown by light pink (C) in Fig. 5, where the overall population remains bounded within [0, 1].We mark this particular kind of oscillations in sub-figure (c), and and in Fig. 5, pink region (C) is shown.The punishers and defectors oscillate periodically, allowing each to dominate the other depending on the time window.Even we are able to detect a very narrow yellow region (K) (See the sub-figures' set (k) in Fig. 5) in the parameter space, where the dynamics of y and z variables exhibit periodic oscillations with period-2.In this region the overall population density x + y + z remains within 1.Note that although society lacks the cooperators, punishers are also special kind of cooperators with some additional power to punish the defectors using their own resources.The absence of cooperators in the community is thus somehow controlled.More importantly, neither the punishers nor the defectors emerge as a dominant winner in this parameter space.Either of them can enjoy being dominant for a time window.However, their fate alters, and the other strategy dominates society for a different time course.Increasing the delay parameter τ leads to an amplification of oscil- lation amplitudes, which undergo period-doubling bifurcation as explained earlier.This results in oscillations with a period of 4 (timeseries for this oscillations are marked in sub-figure (g)), as highlighted in the orange section (G) of Fig. 5.Further increment of τ results in the appearance of eight periodic dynamics (presented in the set of sub-figures (f)) in the deep blue region (F) of Fig. 5, where the proportion of punishers and defectors oscillate.However, these orange and deep blue regions occur only in overcrowded solutions, where the sum of x + y + z exceeds 1.As we continue to increase the value of τ , our model (6) may exhibit chaotic behavior, as indicated by the violet area (E) in the β − τ parameter space.The chaotic phenomenon among the two strate- gies are portrayed by the set of sub-figures (e) in Fig. 5.In this region, the sum of x + y + z may or may not be bounded within 1.Our findings are supported by power spectrum analysis, but these results are not presented here to avoid redundancy.
Figure 5 also displays a beige colored region (D) where the system (6) converges to a steady state with no cooperators.The nature in the player population frequencies for the steady states are shown in the sub-figure (d) in Fig. 5. Furthermore, the system (6) may converge to the extinction equilibrium (0, 0, 0), shown by deep green in Fig. 5, resulting in a population-free state (time series not shown).In the next section, we will perform an analytical calculation for each of these steady states.The white region in the parameter space is identified as a region where the solution may become unbounded or converge to the steady state 0, 0, 1 − ξ σ 3 .This disappearance of attractors is similar to what is discussed in Ref. 46 .In the future, we aim to investigate the underlying cause of this disappearance of attractors in our model.It is worth noting from an ecological perspective that the steady-state solution is now bounded within the range of [0, 1].Interestingly, the cooperator's survivability cannot be facilitated by either the temptation or delay parameter in the entire domain of analysis.Increasing the value of β only helps to promote the density of defectors depending on the initial conditions and other parameter values.The delay parameter leads to the emergence of periodic behavior, allowing the populations to dominate each other regularly.However, periodic and chaotic attractors are not the only possible outcomes that can be obtained by varying the parameters.In Fig. 5, we identify a light purple region where the density of punishers and defectors display quasi periodic oscillations that remain bounded within the range of [0, 1].We additionally identify a regime, shown by cyan (H) in Fig. 5, where both y and z variables oscillate with six periodic dynamics .We present the set of sub-figures (h) to describe the six-periodic oscillation among the competing population Both the deep grey (J) and narrow yellow regions (K) correspond to period-2 oscillation in both the variables, however the deep grey region (J) indicates the overcrowded solution.Time series depictions are identified by sub-figures' set (k), and (j) for the crowded and overcrowded two-periodic oscillatory states respectively.For better visibility, we do mark separate diagrams to portray two different oscillations of the player populations.Period-4 and period-8 oscillations emerge in the orange (G) and deep blue regions (F) , respectively, however again they correspond to overcrowded solution.We represent the population dynamics in order to showcase the four and eight periodic oscillations among the punishers and defectors by sub-figures' set (g), and (f) respectively.Chaotic behavior is identified in the violet area (E) where the dynamics may or may not be bounded within 1. Chaotic behaviors among the two strategies are being shown in the set of sub-figures (e).Cyan (H) and light grey (I) colored regions exhibit period-6 and quasi periodic solutions in both the y and z variables, respectively.Period 6 oscillations are being shown in the set of sub-figures (h), whereas time series showing quasi-periodic oscillations are described in the sub-figures' set (i). Deep green region corresponds to the extinction equilibrium (time series not shown).In the white region either the solution becomes unbounded or the sole existence of defectors is noticed.
traits.Moreover, before the note of a period doubling phenomenon from periodic oscillation to bi-periodic oscillation, we also observe quasi-periodicity among punisher and defector traits, which has been marked by light grey region (I) in Fig. 5. we present the quasi-periodic oscillations in the sub-figures' set (i) of Fig. 5.
In the upcoming subsection, we will shed light on the impact of additional system parameters on our model.Specifically, we will investigate the effects by varying two different parameters while keeping the others fixed.By doing so, we aim to uncover valuable insights into the influence of these specific parameters on our model's behavior.

System parameters and eco-evolutionary dynamics: unveiling the hidden connections.
In the context of our study, it becomes evident that when free space presents a more favorable opportunity for defectors, while keeping other parameter values constant, their survival becomes highly likely.In such scenarios, if the benefit induced by free space towards defectors, denoted as σ 3 , reaches a sufficiently large value, the other strategies face intense competition and may eventually disappear from the societal dynamics.This observation is depicted in Fig. 6a, where we observe that an increase in the mortality rate ξ leads to an equilibrium of extinction, rendering no viable survivors.Conversely, when the mortality rate is moderate, the system exhibits the potential to sustain a society solely composed of cooperators, contingent upon other parameter considerations.Further reducing the mortality rate enables the coexistence of cooperators and defectors.However, an amplification of the free space-induced benefit toward defectors favors their dominance, ultimately eliminating cooperators from the competition.Consequently, under such circumstances, defectors emerge as the sole surviving strategy.
An analogous observation is evident when examining Fig. 6b, which further accentuates this phenomenon of interest.Once again, we observe that as mortality rate escalates, no strategies are able to evolve, ultimately Figure 6.Exploring the role of system parameters in eco-evolutionary dynamics: Subfigure (a) shows that as the mortality rate increases, all strategies face extinction unless free space benefits defectors substantially, leading to their dominance.A significant portion of the parameter space allows for the coexistence of cooperators and defectors, while another portion enables the survival of cooperators alone.Similarly, in subfigure (b), escalating mortality rates result in the extinction of all strategies, except in a region where only cooperators survive.However, this region diminishes with higher temptation.This parameter space exhibits the coexistence of cooperators and defectors, also.In subfigure (c), increasing benefits for cooperators, while keeping other parameters fixed, leads to a surge in their density.Higher mortality rates pose a challenge, but cooperators can persist due to their enhanced benefits.A small region, shown in white, exhibits unbounded dynamics in subfigures (a) and (c).Subfigure (d) illustrates the delicate balance between benefits for defectors and cooperators.When defectors have a greater advantage and exceed the mortality rate, they become the sole survivors.Cooperators can coexist with defectors or thrive as the sole strategy when their benefits surpass both the mortality rate and defectors' benefits.Punishers are absent in the whole figure due to chosen parameters and initial conditions, highlighting the intriguing dynamics.Parameter values: σ 1 = 1.2 , σ 2 = 1.5 , σ 3 = 1.4 , ξ = 1.1 , β = 1.5 , δ = 0.5 , and τ = 0.2 , unless those parameters are varied.Initial conditions for delayed variables: (0.25, 0.25, 0.25), and for non-delayed variables: (0.3, 0.3, 0.3).leading to their eventual extinction over extended time period.However, within the expansive β − ξ parameter space, a distinct region emerges where the survival of solely cooperators becomes feasible.It is important to note, however, that this region diminishes in size as the temptation parameter β increases.This correlation arises from the fact that defectors receive an augmented advantage with a higher β , thereby reducing the space where the sole coexistence of cooperators can only occur.In fact, it is within the moderate range of mortality rate that the majority of the parameter space allows for the coexistence of defectors and cooperators.
In a similar vein, if we delve into the intricate interplay of free space-induced benefits towards cooperators, specifically σ 1 , while keeping σ 3 and other parameters fixed, we anticipate a noticeable surge in the density of cooperators, denoted by x, with each incremental value of σ 1 .The exquisite Fig. 6c illuminates the captivating dynamics that unfold within the system.For lower values of σ 1 , the system may initially lack any cooperators.However, as we venture beyond a critical threshold of σ 1 , contingent upon the values of other parameters, coop- erators reveal their resilience, persisting in the system over extended temporal horizons.While a higher mortality rate poses a formidable barrier to the survival of any individual, the enhanced free space-induced benefits bestowed upon cooperators hold the potential for their sustenance, even under such adversities.Moreover, within the expanse of the σ 1 − ξ parameter space, a significant portion emerges where only defectors can endure due to the relatively diminished contribution of free space towards cooperators.Nevertheless, as the value of σ 1 rises, cooperators resurge, reentering the competitive landscape and forging a coexistence alongside defectors.However, it is crucial to acknowledge the existence of a small region within this parameter space where the dynamics become unbounded over sufficiently long time frames.The intriguing phenomenon of unboundedness manifests itself in a similar fashion within Fig. 6a as well.In the near future, we aspire to delve deeper into this enthralling phenomenon, unveiling the precise mechanisms underlying this unboundedness.
Figure 6d vividly illustrates the intricate interplay between the parameters σ 1 and σ 3 .Notably, when both of these parameters are lesser than the mortality rate ξ = 1.10 , the system reaches an equilibrium of extinction, effectively stabilizing the dynamics in an unwanted scenario.When free space provides defectors with a more favorable opportunity compared to cooperators (indicated by σ 3 > σ 1 ) and σ 3 surpasses the mortality rate ξ , defectors gain substantial advantages, thereby emerging as the sole surviving strategy within that specific parameter space.On the other hand, if σ 1 exceeds both ξ and σ 3 , two distinct scenarios unfold, favoring the survivability of cooperators.In one scenario, the parameter choices allow for the coexistence of cooperators and defectors, contingent upon the specific values of σ 3 .In the second scenario, cooperators thrive as the sole surviving strategy in the societal dynamics, while all others face extinction.It is worth noting that, due to our chosen parameter values and initial conditions, Fig. 6 remains devoid of any punishers, underscoring the fascinating dynamics at play within the system.
Intriguing insights into the behavior of our proposed model await us as we explore the captivating Fig. 7. Remarkably, we discover that the dynamics remain almost unaffected by variations in the delay parameter, despite we have found a profound effect of τ in our study through several earlier discussions (e.g., see Fig. 5).Within Fig. 7a, we embark on a journey to explore the system's response when simultaneously varying two parameters, namely δ and τ , while maintaining fixed values for other parameters at σ 1 = 1.2 , σ 2 = 1.5 , σ 3 = 1.4 , ξ = 1.1 , and β = 1.5 .A discernible pattern emerges from our investigation, shedding light on the interplay between fine imposition and the survival prospects of different strategies.At minimal fine levels, coexistence between punishers and defectors prevails within the societal landscape.However, as we escalate the fine magnitude, defectors too succumb to their inability to thrive, leaving behind a society solely populated by punishers.It is worth noting that these outcomes may exhibit variations when other parameter values are chosen differently.Our selection of parameter values deliberately maintains a moderate temptation to defect ( β = 1.5 ) while endowing punish- ers with a more advantageous position in the free space hierarchy ( σ 2 > σ 3 > σ 1 ).Consequently, cooperators, receiving the least support from free space, struggle to sustain themselves and ultimately face extinction in the long run.Furthermore, this subfigure, Fig. 7a, stands devoid of cooperators, underscoring the inhospitable environment created by the combined influence of fine imposition and delay on their survival prospects.The aforementioned exploration illuminates the intricate dynamics at play, deepening our understanding of the complex interactions within the proposed model.
The captivating Fig. 7b unravels the intricate behavior of our model as we simultaneously vary the parameters ξ and τ , while holding the remaining parameters constant at the same values as in Fig. 7a with δ = 0.5 .Within this visual exploration, we witness the system's remarkable response to these parameter variations, elucidating the delicate balance that governs the survival and coexistence of different strategies.As the mortality rate ξ escalates towards high values, approaching 2, a disheartening scenario unfolds where no strategies manage to endure in the long run.This highlights the adverse consequences of an excessively high mortality rate within the system.However, by reducing the strength of the mortality rate ξ , a fascinating transformation occurs, leading to the emergence of a society comprised solely of punishers.This finding aligns with the insights gleaned from Fig. 7a, highlighting the intricate relationship between mortality rate and strategy dynamics.It is noteworthy to mention that in our proposed model, punishers themselves embody a distinct form of cooperation.Intriguingly, within the two-dimensional parameter space with ξ = 1.0 , we discover a region where cooperators and punishers coexist harmoniously, while defectors remain absent from the societal landscape.This coexistence phenomenon further emphasizes the complex interplay between the investigated parameters.However, as we examine the intricacies portrayed in Fig. 7b, we observe that the delay parameter τ exhibits limited influence on the emerging state within the investigated range.Notably, a thin white portion in the parameter space for ξ = 1.0 unveils the emergence of unbounded dynamics, warranting further investigation into the precise mechanisms behind this intriguing phenomenon.By decreasing the value of the mortality rate ξ even further, we uncover a substantial portion within the parameter space that facilitates the coexistence of all three strategies, fostering biodiversity within the system.As the mortality rate ξ continues to decrease, we witness a captivating transformation where punishers struggle to survive, ultimately leading to a society where cooperators and defectors coexist in delicate harmony.This profound exploration offers invaluable insights into the complex dynamics of our model, unraveling the delicate relationships between mortality rate, delay parameter, and the survival prospects of different strategies.
Both Figs. 6 and 7 showcase numerical simulations based on a fixed initial condition.Specifically, the nondelayed variables are initialized with values of (0.3, 0.3, 0.3), while the delayed variables begin with values of (0.25, 0.25, 0.25).These initial conditions serve as the starting point for investigating the intriguing interplay of parameters in our proposed model and their influence on the emergence of eco-evolutionary dynamics.In the subsequent subsection, we will delve into the pivotal role of initial conditions in our model, analyzing how they shape and contribute to the complex dynamics observed within the system.By examining the impact of different initial conditions, we aim to gain a deeper understanding of the intricate relationships between system parameters, initial states, and the resulting eco-evolutionary dynamics.
Uncovering multistability: insights into eco-evolutionary dynamics.We maintain the parameter values at σ 1 = 1.2 , σ 2 = 1.5 , σ 3 = 1.4 , ξ = 1.1 , β = 1.5 , δ = 0.5 , and τ = 0.2 , consistent with our earlier investigations in Figs. 6 and 7.However, in this subsection, we focus on the role of initial conditions in our model.To explore the range of possible outcomes, we vary the initial conditions (x 0 , y 0 , z 0 ) while satisfying the constraint x 0 + y 0 + z 0 = 0.9 .This constraint allows us to interpret the results from a biological perspective, as x 0 + y 0 + z 0 falls within the range of [0, 1].Similar to our previous investigations, we set the initial condition for the non-delayed variables as (0.25, 0.25, 0.25).Remarkably, we observe that the system converges to different steady states depending on the chosen initial conditions.This finding further reinforces our earlier claim that the emergent dynamics critically depend on the initial conditions.The coexistence of multiple stable states in our model offers valuable insights into the system's stability, robustness, and the possibility of regime shifts.In particular, in Fig. 8, we present the results that reveal the existence of four distinct attractors.These attractors represent different stable strategies that can persist over time.Under specific conditions, we can observe the coexistence of punishers and defectors within the system.Moreover, depending on the initial conditions, the system may exhibit stable states consisting solely of cooperators, punishers, or defectors.The identification of these multiple stable states sheds light on the complex dynamics and potential outcomes in our eco-evolutionary model.Understanding the interplay between initial conditions and the resulting steady states is crucial for comprehending the long-term behavior and potential shifts in ecological and evolutionary systems.
In the upcoming section, we derive the steady states and analytically determine the point of Hopf bifurcation in our model.Through rigorous analysis and numerical validation, we establish the accuracy of our findings.

Analytical results
Various steady states and their biological relevance.The proposed model ( 6) results in eight steady states.They are briefly discussed in the following, • The extinction state S 0 : The state of extinction, denoted as S 0 , corresponds to the equilibrium point (0, 0, 0) in the dynamical system.At this state, each of the player populations eventually die out due to intense competition over long periods of time.• Punisher and Defector-free state S 1 : In this steady state, only individuals with cooperative strategy have the opportunity to survive, leading to the extinction of all other types of players in the game.The stationary point associated with this state is given by σ 1 − ξ σ 1 − 1 , 0, 0 .Under these circumstances, cooperators have the ability to survive and ultimately dominate the entire game.• Cooperator and defector-free state S 2 : This state exclusively enables the survival of punishers, resulting in the eventual extinction of all other populations in the game.The stationary point associated with this state is 0, σ 2 − ξ σ 2 − 1 , 0 .In this scenario, punishers gain the most significant advantage and can solely dominate the entire game.However, in an interactive context, if a society consists of exclusively punishers, they would act as cooperators because there would be no defectors to punish.• Cooperator and punisher-free state S 3 : This state allows only defectors to survive in a society, which supports the fundamental theory of the prisoner's dilemma game, as defectors receive the greatest benefit compared to all other population strategies.The equilibrium point associated with this state is 0, 0, 1 − ξ σ 3 .In such scenario, where only one population strategy can survive in the long run, there will be no interactive action between different player populations with different strategies.• Defector-free state S 4 : In this steady state, cooperators and punishers get the chance to survive side by side, with no defectors in the interaction.The equilibrium point gets the form (x * , y * , 0) , where, Here, the defectors are unable to survive, whereas, the cooperators and the punishers jointly survive.Punishers also behave like cooperators in such state.• Punisher-free state S 5 : This state in the proposed model represents the basic scenario of the prisoner's dilemma game.In the long run, only the cooperators and defectors interact with each other, with no punishers present in the society.In such a scenario, the chances of implementing cooperation in different ways are reduced because the defector is not constrained in dealing with cooperators in the absence of punishers.The stationary point can be expressed as (x * , 0, z * ) , where • Cooperator-free state S 6 : This steady state operates similarly to the previous state, S 5 , where, in the long run of the dilemma, the punishers and the defectors can interact with each other.However, unlike state S 5 , the defectors in this state face a slight constraint when interacting with the punishers.There would be a slight reduction in their payoffs when interacting with punishers than the cooperators.On the other hand, the punishers also behave like cooperators, but when interacting with the defectors, they face a loss of the fine δ to punish t he defe c tors.The e qui libr ium p oint is g iven by (0, y * , z * ) , w here , and .

Figure 8.
Exploring the impact of initial conditions on eco-evolutionary dynamics: Varying initial conditions within the constraint x 0 + y 0 + z 0 = 0.9 , we uncover diverse steady states, highlighting the profound influence of initial conditions on emergent dynamics.This multistability provides valuable insights into system stability, robustness, and regime shifts.Notably, our findings reveal four distinct attractors representing different stable strategies, including coexistence of punishers and defectors.Depending on the chosen initial conditions, the system can exhibit stable states with cooperators, punishers, or defectors alone.These discoveries deepen our comprehension of the intricate dynamics and potential outcomes in our eco-evolutionary model, underscoring the critical interplay between initial conditions and steady states.Parameter values: σ 1 = 1.2 , σ 2 = 1.5 , σ 3 = 1.4 , ξ = 1.1 , β = 1.5 , δ = 0.5 , and τ = 0.2 .Initial conditions for delayed variables: (0.25, 0.25, 0.25).
• State of co-existence S 7 : This steady state allows all three types of player populations to coexist and interact with each other simultaneously.It is considered the most biologically significant among all eight states as no player population goes extinct in the long run.The equilibrium point for this state is (x * , y * , z * ) , where Tracking the point of occurrence of Hopf bifurcation.Assuming (x * , y * , z * ) as the coordinates of a specific steady state, we explore the occurrence of the Hopf bifurcation by progressively increasing the time delay variable τ imposed on the system from this stable steady state.We refer to the value of τ at which the steady states begin to lose stability as the critical value of τ or τ c .The linearization form of the proposed system ( 6) is presented below.
By calculating the linearized version of the above system about the general steady state (x * , y * , z * ) , we obtain where The characteristic equations of the linearized system (7) are obtained as Ṡ = JX , where In this case, J represents the Jacobian of the linearized system of equations.To obtain its characteristic equation, we set det(J − I) = 0 , where is an eigenvalue of the Jacobian matrix J. From the resulting eigenvalues, we can deduce the critical value of τ .Therefore, we have the following transcendental equation from det(J − I) = 0: where To simplify the transcendental equation, we make the following substitution: (τ ) = a(τ ) + ib(τ ) and (τ c ) = a(τ c ) + ib(τ c ) .We then separate the real and imaginary parts of the equation by using the transforma- tion e a±ib = e a (cos(b) ± isin(b)) on Eq. ( 9).This yields � .
( Using the relation in Eq. ( 11), we can determine the value of the delay time variable τ c at which the Hopf bifurcation occurs.At this value of τ , the eigenvalues of the Jacobian become purely imaginary.Therefore, we set a(τ c ) = 0 .This provides the following form of Eq. ( 11): Let us consider that Consequently, we find the following relations, Using these two relations, we can easily conclude Since the left-hand side of Eq. ( 13) includes the terms (A − B) 2 and (C − D) 2 , and considering the assump- tions we have made about these values, we can substitute them and rewrite Eq. ( 13) as follows: After performing extensive calculations, we can derive the following relationship by multiplying cos(2b * τ c ) with the first equation of system (12) and sin(2b * τ c ) with the second equation of the same system.
( 16) brings forth one of the solutions as cos(b * τ c ) = 0 .Therefore, we need to determine the value of b * for which the relation cos(b * τ c ) = 0 holds true.We can use this relationship to calculate the critical value τ c , at which Hopf bifurcation occurs.Therefore, the value of τ obtained from this relationship can be considered as the value of τ c .By modifying Eq. ( 14), we can obtain a sixth-degree polynomial in terms of b * .Substituting the value of sin(b * τ c )| cos(b * τ c )=0 from Eq. ( 15), we can derive an equation in terms of b * , which is sufficient to calculate τ c .The resulting sixth-degree equation is where The roots of Eq. ( 17) are the six values of b * such that the relation for the critical value of τ holds, i.e., cos(b . For n = 0 , the critical value τ c can be calculated as τ c = π 2b * , where the value of the required b * can be calculated from Eq. ( 17).Note that one can obtain at most six different real values for b * .Only one of these real roots provides us with our desired critical value of the delay parameter at the Hopf bifurcation.
Validation of analytical finding.In Fig. 3, we observe the occurrence of Hopf bifurcation, where a stable equilibrium point becomes unstable, giving rise to a periodic solution at τ = 8.98 .This critical point is consist- ent with the analytically derived τ c .Furthermore, in Fig. 9, we demonstrate the bifurcation diagram of our system (6) by varying the delay parameter τ within the range [2,3], using a fixed step length of 0.002, with ξ = 0.40 , β = 1.50 , δ = 0.56 , σ 1 = 1.20 , σ 2 = 1.70 , and σ 3 = 1.0 .Despite the fact that free space favors punishers due to σ 2 > σ 1 > σ 3 , punishers eventu- ally become extinct in the long run (as seen in Fig. 9c).In the steady-state regime, defectors can dominate cooperators, even though they are less favored by the free space (as seen in Figs.9b,d).However, this steady-state loses its stability at τ = 2.706 , resulting in a limit cycle oscillation.The occurrence of Hopf bifurcation at this point in the system parameter space confirms our analytically derived τ c = 2.706 .The maximum Lyapunov exponent 1 , shown with solid purple line in subfigure (a), remains at zero beyond this bifurcation point, confirming the periodic oscillation of the system in the range studied.Additionally, the second-largest Lyapunov exponent 2 , shown with cyan dashed line, touches zero at this bifurcation point, further validating our bifurcation diagram.However, beyond a certain value of the delay parameter τ , the system's solution becomes overcrowded, as the values of (x + y + z) surpass the upper bound of 1.We highlight this ecologically feasible range by plotting a horizontal dash-dotted line in subfigures (a) and (d), and solutions below this line remain ecologically meaningful and interpretable.

Discussion
Eco-evolutionary dynamics investigates how ecological and evolutionary processes interact and influence each other over time.It involves understanding the feedback loops between ecological and evolutionary processes and how they shape the dynamics of populations and communities.Ecological processes such as competition, predation, and resource availability can drive natural selection and thus influence the evolution of traits in populations.Conversely, the traits and genetic diversity within populations can shape ecological processes and affect the structure and functioning of communities.Eco-evolutionary dynamics can have important implications for conservation biology, as changes in ecological processes or environmental conditions can alter the direction or rate of evolution and vice versa.For example, human-induced environmental changes such as habitat fragmentation or climate change can alter the selective pressures on populations, leading to changes in their evolutionary trajectory.Overall, understanding the complex interplay between ecological and evolutionary processes is crucial for understanding the functioning and resilience of ecosystems in the face of environmental change.
Our present study examines the interactions and long-term effects between ecological and evolutionary processes by considering a classic game, viz. the prisoner's dilemma.The prisoner's dilemma is often used as a metaphor for situations in which two individuals or groups must choose between cooperation and competition and in which self-interest can lead to a suboptimal outcome for both parties.It has important applications in economics, political science, and evolutionary biology.Our study suggests free space being an ecological variable, can significantly influence the prisoner's dilemma game's evolutionary dynamics.Free space here acts as an individual that benefits another individual at a cost.They lose their fitness by helping others without any expectation of reward.Although this form of philanthropic activity may seem to be foolish and may raise a question on the evolutionary basis of altruism among researchers, one can still observe such altruistic behavior in a variety of forms, such as donating to charity, volunteering time to help others, or even sacrificing one's own life to save another.Altruism can also occur in non-human species, such as when a mother animal risks her safety to protect her offspring.We incorporate this vital aspect of social behavior which has important implications for fields such as psychology, sociology, and evolutionary biology.In contrast, punishment is an effective strategy to discourage individuals from engaging in undesirable behavior by making the cost of that behavior exceed the benefit.Punishment can provide a sense of justice or closure for those harmed by the conduct and can prevent people or groups from behaving in ways that are undesirable or unacceptable.To incorporate punishment into our game, we introduce a fraction of punishers.This approach, combined with the classical prisoner's dilemma game and two additional strategies, helps us to understand how ecological and evolutionary processes interact and impact population and community dynamics.However, despite making progress, we have limited knowledge of how time lag affects eco-evolutionary game dynamics.The study of dynamical systems with delay is an active research area in mathematics, physics, and engineering.It has practical applications in fields such as control theory, signal processing, and neuroscience.Through numerical simulations, we find that delay can lead to oscillation in our system, which is impossible in non-delayed cases with those sets of parameter values.The evolution of strategies depends on several factors, shown with solid purple line becomes zero at the bifurcation point and remains unchanged beyond the critical point.The second largest Lyapunov exponent 2 (cyan dashed line) touches zero at the bifurcation point and remains negative for other values of τ , and the third exponent 3 (pink dotted line) remains negative throughout the range of τ .Beyond a certain value of τ the system exhibits overcrowded solution as the total population x + y + z exceeds the maximum value 1.A green horizontal dash dotted line is plotted to mark the boundary of the overcrowded solution from that of the bounded solution.(b-d) Showcase the effect of delay τ on each of the population variables.Punishers are unable to survive for the chosen parameter values even though they receive most benefit from the free space compared to the cooperators and defectors.We choose initial fractional quantities for the non-delayed variables to be (0.1, 0.2, 0.5) and for the delayed variables to be (0.3, 0.3, 0.3).such as the resources punishers use to fine defectors and the charitable contribution of free space.The survival of players using different strategies also depends on the initial fraction of the population.The system can settle into various stable states, and understanding multistability can provide valuable insights into the complex behavior of our system.We derive all steady states and calculate the point of Hopf bifurcation, where the steady state loses its stability, and delay causes the system to undergo periodic oscillations.This periodic attractor enables each strategy to dominate the other in a cyclical manner under suitable circumstances.Although defection appears to be the rational choice in the classical prisoner's dilemma game, the intricate interactions between ecological and evolutionary processes prevent any strategy from dominating the others in the long run.Our proposed eco-evolutionary model can behave like a three-species cyclic dominance system for an appropriate selection of parameters, with each species dominating and being outcompeted by the next.Following the Hopf bifurcation, the system may undergo a series of period-doubling bifurcations as the parameter value increases, resulting in a successive doubling of the period of oscillation and the emergence of chaotic behavior.However, this can lead to overcrowding and limit the biological implications of our study.Our study also has several limitations worth exploring further.Our study uses paired communications to characterize group interactions and is conducted for the prisoner's dilemma game.Investigating the impact of higher-order interactions 81,82 and how the outcome varies for other evolutionary games 83,84 might be fascinating.In conclusion, our analysis provides valuable insights into the complex eco-evolutionary dynamics and contributes to a better understanding of group decision-making and the emergence of moral behavior in multidimensional social systems.We hope that our study will inspire further research in this area and facilitate the development of effective strategies for managing and conserving ecological communities. https://doi.org/10.1038/s41598-023-41519-1 https://doi.org/10.1038/s41598-023-41519-1

Figure 3 .
Figure 3. Bifurcation diagram and the Lyapunov exponents of the proposed delay model: (a)The effect of the delay parameter τ on the dynamics of the total population x + y + z of all three species, that indicates the emergence of periodic oscillation via Hopf bifurcation at τ = 8.99 which then leads to period doubling bifurcation as delay increases and eventually to chaos.The three Lyapunov exponents 1,2,3 of the delayed system are also displayed among which 1 (solid purple) is the maximum that validates whether the dynamics is a steady state, periodic or chaotic oscillation contingent to its value being negative, zero or positive, respectively.The second largest exponent 2 (cyan dashed) remains negative throughout the delay interval except at the bifurcation points where it takes zero value and 3 (pink dotted) remains negative throughout the interval.The individual bifurcation diagrams are also shown in (b) for cooperators x, (c) for punishers y, and (d) for defectors z, respectively.The following parameter values are chosen: ξ = 1.20, β = 1.60, δ = 0.30, σ 1 = 1.35, σ 2 = 1.50, σ 3 = 1.35 and the initial conditions are set at (0.1, 0.2, 0.5) for the non-delayed variables and (0.3, 0.3, 0.3) for the delyaed variables.The figures are plotted by taking the τ step length of 0.005 for the bifurcation diagrams and 0.04 for the Lyapunov exponents in the range[8,12].From the individual bifurcation diagrams it is clear that the chosen parameter values correspond to the coexistence of only the punishers and defectors, with punishers being the dominating population.

Figure 4 .
Figure 4. Different dynamical structures of the player populations along with their power spectra: The effect of the temptation parameter and delay parameter on the temporal dynamics of the delayed system are displayed for three different choices of (β, τ ) values in three different panels.Time series of the cooperators, punishers and defectors are portrayed in the first, second and third columns, respectively.The fourth column corresponds to the trajectory in the y − z parameter space, and finally, in last column the power spectra using fast Fourier transformation is presented for three different temporal dynamics.First panel (a-e) exhibits the chaotic dynamical behavior in y and z variables for the chosen parameter values (β, τ ) = (1.93,11.64) , second panel (f- j) corresponds to the emergence of two periodic oscillation in punishers' and defectors' dynamics for the chosen parameter values (β, τ ) = (1.64,10.73) , and third panel (k-o) displays the quasi-periodic dynamics among the punisher and defector populations with irregular periodicity for the parameter values (β, τ ) = (1.90, 10.50) .In all three cases, the cooperators remain extinct due to the choice of the other parameter values which are fixed at ξ = 1.2, δ = 0.3, σ 1 = 1.35, σ 2 = 1.5, σ 3 = 1.35 .These three distinct dynamical states are also justified by the power spectrum analysis that unveils the nature of the peaks of the power spectra measured with respect to the normalized frequency.The initial values of the population fractions are taken to be (0.1, 0.2, 0.5) for the nondelayed variables, and (0.3, 0.3, 0.3) for the delayed variables.We have observed the resultant dynamics after wiping out the initial 98 × 10 5 out of the 10 7 numbers of iterations.

Figure 5 .
Figure5.Two dimensional parameter space due to the simultaneous interplay of the temptation parameter β and the delay parameter τ : The β − τ parameter space unveils various temporal dynamics of different population variables in the region β ∈ [1, 2] and τ ∈[8,12] .The other parameter values are fixed at ξ = 1.20, δ = 0.30, σ 1 = 1.35, σ 2 = 1.50, and σ 3 = 1.35, the initial conditions are set to be (0.1, 0.2, 0.5) and (0.3, 0.3, 0.3) for the non-delayed and delayed variables, respectively.Due to this particular choice of the parameter values both the cooperators and the defectors receive equal benefit from the free space (i.e., σ 1 = σ 3 ), however we end up with a society free from cooperators in the parameter region under study because of the initial advantage to the defectors over the cooperators.Time series of the punishers and defectors corresponding to different colored regions in the parameter space are shown, cooperators are not shown since they remain extinct throughout the space.The mauve region (A) on the left corresponds to a community only consisting of punishers that oscillate periodically.Time series in order to show oscillations in the player frequencies for this region are depicted by sub-figure (a).Periodic coexistence of both the punishers and defectors are observed in both the regions marked with lime yellow (B) and light pink (C) , but the region (B) corresponds to the overcrowded solution as y + z > 1 , and the region (C) corresponds to the bounded solution.Sub-figure (c) highlights the periodic oscillatory time series data depicted by the of the region (C), whereas sub-figure (b) showcases the overcrowded scenario, marked by (B) region.Steady state coexistence of both punishers and defectors is noticed in the beige region (D) .The steady state time series data is marked in the sub-figure (d).Both the deep grey (J) and narrow yellow regions (K) correspond to period-2 oscillation in both the variables, however the deep grey region (J) indicates the overcrowded solution.Time series depictions are identified by sub-figures' set (k), and (j) for the crowded and overcrowded two-periodic oscillatory states respectively.For better visibility, we do mark separate diagrams to portray two different oscillations of the player populations.Period-4 and period-8 oscillations emerge in the orange (G) and deep blue regions (F) , respectively, however again they correspond to overcrowded solution.We represent the population dynamics in order to showcase the four and eight periodic oscillations among the punishers and defectors by sub-figures' set (g), and (f) respectively.Chaotic behavior is identified in the violet area (E) where the dynamics may or may not be bounded within 1. Chaotic behaviors among the two strategies are being shown in the set of sub-figures (e).Cyan (H) and light grey (I) colored regions exhibit period-6 and quasi periodic solutions in both the y and z variables, respectively.Period 6 oscillations are being shown in the set of sub-figures (h), whereas time series showing quasi-periodic oscillations are described in the sub-figures' set (i). Deep green region corresponds to the extinction equilibrium (time series not shown).In the white region either the solution becomes unbounded or the sole existence of defectors is noticed.

Figure 7 .
Figure 7. Unveiling the interplay of delay, mortality rate, and punishment in curbing selfish behavior of defectors: (a) Corresponds to the exploration of the response to simultaneous variations of δ and τ , uncovering a pattern of strategy survival shaped by fine imposition.Coexistence between punishers and defectors prevails at minimal punishment while escalating punishment leads to the dominance of punishers and the extinction of defectors.These outcomes depend on our specific parameter values, emphasizing the delicate balance between temptation, hierarchy, and strategy survival.(b) Unveils the intricate behavior of our model when varying ξ and τ , shedding light on the interplay between mortality rate and strategy dynamics.High mortality erodes all strategies, but reducing ξ allows for a society of sole punishers.Coexistence between cooperators and punishers emerges in a very narrow parameter region, free from defectors.The influence of the delay parameter τ is limited within the investigated range, except for a thin white portion with unbounded dynamics.Decreasing ξ further leads to the coexistence of all strategies, fostering biodiversity, and eventually, to a society of cooperators and defectors.Parameter values: σ 1 = 1.2 , σ 2 = 1.5 , σ 3 = 1.4 , ξ = 1.1 , β = 1.5 , δ = 0.5 , and τ = 0.2 , unless those parameters are varied.Initial conditions for delayed variables: (0.25, 0.25, 0.25), and for non-delayed variables: (0.3, 0.3, 0.3).

Figure 9 .
Figure 9. Verification of the analytically obtained value of the critical time delay τ c : (a) Demonstrates the bifurcation diagram and three Lyapunov exponents of the delay model by varying the delay parameter τ in the range[2,3].We use the following parameter values: ξ = 0.40 , β = 1.50 , δ = 0.56 , σ 1 = 1.20 , σ 2 = 1.70 , and σ 3 = 1.0 .For the chosen parameter values the system undergoes a Hopf bifurcation at the critical point τ c = 2.706 , that matches perfectly with the analytically derived τ critical value.The largest Lyapunov exponent 1 shown with solid purple line becomes zero at the bifurcation point and remains unchanged beyond the critical point.The second largest Lyapunov exponent 2 (cyan dashed line) touches zero at the bifurcation point and remains negative for other values of τ , and the third exponent 3 (pink dotted line) remains negative throughout the range of τ .Beyond a certain value of τ the system exhibits overcrowded solution as the total population x + y + z exceeds the maximum value 1.A green horizontal dash dotted line is plotted to mark the boundary of the overcrowded solution from that of the bounded solution.(b-d) Showcase the effect of delay τ on each of the population variables.Punishers are unable to survive for the chosen parameter values even though they receive most benefit from the free space compared to the cooperators and defectors.We choose initial fractional quantities for the non-delayed variables to be (0.1, 0.2, 0.5) and for the delayed variables to be (0.3, 0.3, 0.3).