Evidence of Chaos in Eco-epidemic Models

We study an eco-epidemic model with two trophic levels in which the dynamics is determined by predator-prey interactions as well as the vulnerability of the predator to a disease. Using the concept of generalized models we show that for certain classes of eco-epidemic models quasiperiodic and chaotic dynamics is generic and likely to occur. This result is based on the existence of bifurcations of higher codimension such as double Hopf bifurcations. We illustrate the emergence of chaotic behavior with one example system.


Introduction
Mathematical models are essential tools in order to understand the mechanisms responsible for persistence or extinction of species in natural systems.In ecological models persistence is in general desired.By contrast, investigations in epidemic models usually aim at finding mechanisms that lead to the extinction of the parasites or infections (e.g.[39]).However, it is known that diseases can not only greatly affect their host populations, but also other species their host populations interact with (e.g.[2]).
Most of the existing models in eco-epidemiology consider a disease in the prey population [5,14,58].There are only a few models where the predator population is infected [60].Some of the latter can exhibit sustained oscillations which are absent from the uninfected ecological model under consideration [2,62].
In ecology as well as in epidemiology oscillations are associated with destabilization.The reason is that extinction of the population due to natural fluctuations becomes very likely when the oscillation drives the population to low abundances [19,44,47,56].Such extinction events are less critical if a species is spatially separated in several subpopulations.The subpopulations allow for a repopulation by migration after the extinction of one subpopulation.However, in such a case synchronous dynamics between the subpopulations due to the coupling by migration can have a devastating effect: it increases the possibility of global extinction of the whole species [31,37] when all subpopulations have a minimum at the same time.In epidemiology such a synchronization of the dynamics can again be desired and induced by pulse vaccinations [22] in order to cause the extinction of the disease.
In some models where the prey is infected also chaotic long-term dynamics have been observed [13].In contrast to oscillations, the effect of chaos on the stability in ecological models is a question of debate for a long time [25,30,50].A popular view has been that chaos has a destabilizing effect because of the associated boom and burst dynamics [11].However, for a population that consists of subpopulations as described above it has been shown that chaotic dynamics can be of advantage in order to stabilize the whole population.Although chaotic dynamics increase the number of local extinctions of subpopulations it reduces the degree of synchrony between different patches.Consequently it reduces the probability of a global extinction [1,49].Thereby it seems essential that the subpopulations are chaotic in isolation [20].With the coupling the subpopulations can show simpler dynamics but the subpopulations tend to be out of phase.In this sense it can be of advantage for the total population if the subpopulations exhibit chaotic dynamics or are close to chaotic parameter regions.This observation may explain why some ecological systems appear to be at the edge of chaos [54].Further, diffusion-induced complex dynamics as been found in continous spatial predatorprey systems [7,45,55].But also isolated populations can exhibit persistent chaotic dynamics.Benincà et al. [10] observed chaotic dynamics of a food web in laboratory mesocosm that last for more than 6 years.From an evolutionary point of view ecological models should evolve into chaotic parameter regions if chaotic dynamics are of advantage for the population persistence.
To study stabilizing or destabilizing effects bifurcation theory is often applied to find the parameter regions where the steady state for the ecological or epidemiological model is stable with respect to perturbations.These analyses are based on specific models describing the relevant processes in form of ordinary differential equations (ODEs).Recently another approach has been developed which is based on generalized models where the exact mathematical form of the processes entering the right hand side of the ODEs is not specified [27].In spite of the fact that the mathematical functions are not known in detail, it is possible to analyze the stability properties of the steady state and draw conclusions about possible destabilization mechanisms of the steady state when a system parameter is varied [26], as well as the emergence of chaotic dynamics [25].
In this paper we use the concept of generalized model to study the impact of a disease of the predator on the dynamics of a predator-prey system.To this end we couple a generalized ecological model with an epidemic one and study the stability properties of the steady state.The advantage of the method is that our results are not restricted to a particular model but apply to certain classes of models since we use parameters encoding the shape of predatorprey functional responses and incidence functions as bifurcation parameters.We find that the coupling of both models introducing a disease in the predator population leads to complex dynamics such as quasiperiodic and chaotic motion.This complex dynamics are solely based on the interplay of demographic and epidemic modeling since the demographic and the epidemic model alone are not capable to exhibit chaotic motion.As a result we show that chaos is generic and prevalent for certain classes of eco-epidemic models.
The paper is organized as follows: In Sec. 2 we discuss the demographic as well as the epidemic model and introduce the eco-epidemic model.Furthermore, we normalize the model and discuss the possible emergence of bifurcations for the steady state.In Sec. 3 we analyze the stability properties of the steady state for the predator-prey model with and without the disease in the predator population and compare the results.As a conclusion we obtain classes of systems in which we expect to find quasiperiodic and chaotic behavior.Since these findings are based on mathematical theorems we can only state the existence of parameter regions where the dynamics is chaotic.To find out the size of these parameter regions and the hence, their relevance for the dynamics of the eco-epidemic model we investigate a specific model in Sec. 4 to show explicitly the emergence of chaotic dynamics.Finally we summarize the results in Sec. 5.

The generalized eco-epidemic model
Before we begin our analysis we will briefly outline the construction of one of the simplest predator-prey models and of one of the most elementary epidemic models, which together form the building blocks for the more general ecoepidemic model we would like to consider.
The basic demographic model in general accounts for two interacting species.The nature of interactions can be of competing, predator-prey or symbiotic nature.A typical formulation for a predator-prey model is given by where S is the specific growth rate of the prey X. Apart from predation the growth of the prey is limited by intraspecific competition assumed to increase quadratically in X with the coefficient M X thereby giving a logistic evolution for the prey dynamics.The predation is expressed by the so-called per capita functional response G(X).The efficiency of biomass conversion is given by the yield constant E. M Y is the mortality rate of the predator population.
This simple model structure has been analyzed for a variety of different functional responses G(X) describing the prey-dependent predation rate (e.g.[36]).However, in our case the function G(X) is not specified in order to keep the model more general.In Sec. 3 we will discuss some properties of this underlying ecological model.
Classical epidemic models partition the population into several epidemiological classes, for a thorough review see [33].The population, in our case the predator population Y , is usually split into susceptibles Y S , infected Y I , and recovered Y R .The latter may be thought of as being immunized, at least for a period of time, after which they return into the class of susceptibles.Following this population division, in absence of vital dynamics, i.e. demographic terms to account for births and natural deaths, a simple SIRS model would be written as assuming linear transition rate γ from the infected to the recovered class.The recovered become susceptible again with a fixed rate δ.

No. Functional form
Comments Citations non-monotone incidence [61] 9 ∼ S I A+S+I asymptotic incidence [17,46] Table 1 Different proposed functional forms for the incidence function In general, pathogen transmission is expressed by interactions among individuals.The latter are modeled by the incidence function λ(Y S , Y I ), for which the most common approaches are the mass action λ(Y S , Y I ) = bY S Y I and the so-called standard incidence function or frequency-dependent transmission λ(Y S , Y I ) = bY S Y I /(Y S + Y I ).In both cases susceptibles Y S and infected Y I are assumed to be well-mixed and hence, to interact randomly.However, it is not clear if the assumption of random interactions and an equal distribution of infected and uninfected is appropriate to describe pathogen transmission in wild populations.Both, small-scale experiments as well as observed disease dynamics, give evidence that simple mass action is not an adequate model in many situations [43].The simplest argument for an asymmetry of the incidence function is that due to a patchiness in the disease on average each infected individual is more likely to have an infected neighbor.The more biological details are taken into account, the more complex the incidence function may be.For instance, Capasso and Serio [12] introduced a saturated incidence function.Such a saturation can be caused by crowding effects at high infection levels or by protection measures the susceptible individuals take.Liu et.al. [40] proposed a more general incidence function of the form λ(Y S , Y I ) = kY S Y p I /(1 + mY p I ).Additionally a variety of other incidence functions have been investigated by various authors (e.g.Table 1).A universal approach has not been found yet.However, using the generalized approach we avoid to specify the incidence function but study more generic properties of the model class under consideration.
To analyze the effect of a disease in the predator population Y on the dynamics of interaction with the prey X we combine the demographic model Eqs.(1) with the SIRS epidemic model Eqs.(2) as follows Here we assume neither vertical transmission, nor vertical immunity, i.e. that the infected predators as well as the recovered predators reproduce only susceptibles.Furthermore, we suppose that the disease can, in principle, influence the demographic parameters.The disease may induce a disease related mortality rate µ and reduce the predation and reproduction rates of the infected Y I expressed by the factors α and β respectively.
All parameters and terms denoted by Greek letters are related to the disease.Note, that in the absence of weakening effects of the disease (i.e.α = β = 1 and µ = 0) the combined model Eq.( 3) reproduces the population dynamics of the uninfected model Eqs.(1) , with Y = Y S + Y I + Y R .This means that the disease could have in principle no influence on the ecological dynamics.
To analyze the dynamics of model ( 3) one would start by computing the steady state and its stability with respect to perturbations.But a local stability analysis cannot be performed since an analytical computation of the steady states is impossible because G(X) and λ(Y S , Y I ) are not specified but assumed to be general functions.However, this difficulty can be overcome using the normalization procedure for the generalized models described in [27].To use this approach we assume that a positive steady state (X * , Y * S , Y * I , Y * R ) exists.
We now define normalized variables x := X/X * , y The details of the normalization and the definitions of the newly introduced scale parameters a i , b, b, f α , fα , f β , fβ , e s , ẽs , m x , mx , m y and my are given in the Appendix.As an advantage of this approach these parameters are easy to interpret in the biological context.The scale parameters a x , a s , a i and a r for instance encode the inverse timescales of the normalized state variables.They measure the relation between the lifetimes of the different species.All other scale parameters are between 0 and 1 and describe weight factors of certain processes of the model at the steady state.
The losses due to intraspecific competition relative to the total losses within the prey are represented by the parameter m x .To be specific, if m x is close to 1 the losses of prey due to intraspecific competition preponderate.The parameter mx = 1 − m x expresses losses caused by predation.f α and fα are the fractions of prey consumed by infected predators and healthy predators respectively.In the same way b is the fraction of healthy predators that are susceptible and b = 1 − b the fraction of healthy predators that are recovered.Further, the parameter e s represents the weight factor of the natural growth terms of susceptibles due to consumption of X.At the steady state the fraction of gains due to recovered predators that become susceptible again is given by ẽs = 1 − e s .The natural mortality for the predator relative to the total losses is expressed by m y .
In the normalized model the steady state under consideration is known (x * , y * ) = (1, 1).The stability of this steady state depends on the eigenvalues of the Jacobian.The steady state is stable if all eigenvalues have a negative real part.Consequently only two bifurcations can separate stable from unstable parameter regions: the saddle-node type bifurcation where a real eigenvalue crosses the imaginary axis and a Hopf bifurcation where a pair of complex conjugate eigenvalues crosses the imaginary axis.
Because all normalized state variables and the normalized processes (l(y s , y i ), g(x)) are equal to one at the steady state, the Jacobian of the normalized model contains in addition to the scale parameters only the derivatives of the normalized processes in the steady state.We define as the generalized parameters.These parameters can be interpreted as nonlinearity measures of the corresponding functions with respect to the variable of the derivative.If the function G(X) is linear in X the derivative of the normalized function g x is equal to one.It is zero for a constant function and two for a quadratic function.To be consistent with previous publications we let g x be the predator sensitivity to prey [24,25,27].In the same sense we denote by l s and l i the incidence sensitivity to susceptibles and to infected respectively.
In summary the Jacobian consists of 10 scaling parameters and 3 generalized parameters.How to obtain the test functions for the above mentioned bifurcations from the Jacobian is described in detail in [26].These test functions enable us to draw 3D bifurcation diagrams as described in detail in [52].
Since we are essentially interested in the influence of different mathematical expressions for the functional response and for the incidence function, we focus our bifurcation analysis on the generalized parameters g x , l s and l i .We chose the other scale parameters according to biological reasoning.It is known that in many cases the timescale for the lifetime of species belonging to different trophic levels slows down with each higher trophic levels [32].Hence, we could assume that the inverse timescale of the susceptible predators is less than half the timescale of the prey, i.e. a s = 0.4a x .By renormalizing the timescale we can say that a x = 1 and a s = 0.4.It is further reasonable to expect that the timescale of the infected predators is slightly larger than the timescale of the susceptible predators since we suppose that their overall lifetime is shorter.Let us assume a i = 0.5.Clearly, this intuitive way is much more appropriate than guessing some abstract parameters.If we would analyze a specific real system at the steady state we could, in principle, also gain an appropriate value for each scale parameter by measuring the corresponding rates.Approximating all other scale parameters, we end up with three parameters which we consider as the most interesting bifurcation parameters.These are the sensitivity of the predator with respect to prey g x and the sensitivity of the incidence function with respect to susceptibles l s and infected l i .The computation of 3D bifurcation diagrams allows us to discuss the stability properties of the eco-epidemic model depending on the mathematical form of the functional response G(X) and the incidence function λ(Y S , Y I ).
3 Stability of the steady state: from local to global bifurcations

Absence of diseases
Before we analyze the effect of an infection on the predator-prey interactions, we take a look at the generalized predator-prey model in the absence of infected individuals.It is known that the predator-prey system Eq.( 1) can exhibit self-sustained oscillations if the functional response G(X) is nonlinear in X for instance a Holling type II function [36].A typical example would be the Rosenzweig-McArthur model [48].These oscillations appear due to a supercritical Hopf bifurcation.Figure 1 shows the bifurcation diagram of the generalized predator-prey model given by Eq. (1).As mentioned above a x is set to be one so that a y corresponds to the relative inverse timescale.We see two saddle-node type bifurcation surfaces (blue) and one Hopf bifurcation surface (red).The steady state is stable in the top volume of the diagram.If one of the bifurcation surfaces is crossed due to a parameter variation the steady state becomes unstable.The only biologically sound parameter range for the scale parameters a y and m x lies between [0,1] since both parameters express some kind of weight factor measured in relation to other scale factors.
Firstly note, that such a destabilization never occurs for a variation of a y .The timescale has therefore no influence on the stability of the steady state in this model class.
Secondly, the Hopf bifurcation surface exceeds the biologically relevant parameter range for g x ≥ 1.For this reason, Hopf bifurcations cannot be found in this model class if g(x) and therefore G(X) are linear functions (g x = 1).This situation corresponds to the Lotka-Volterra model coupled with logistic growth.However, from a biological perspective, models should allow lower values of the sensitivity to prey, i.e. g x < 1. Due to a limited consumption of prey the functional response G(X) should saturate at high amounts of prey.Since saturation is related to low values of g x , Hopf bifurcations should likely occur in biological realistic models.
Even, lower values of the sensitivity to prey, i.e. g x ≤ 0, are rather unlikely and occur only in systems with non monotonic functional responses.Biologically this region where the predation decreases with increasing prey can be related to inhibition effects or group defense techniques of the prey [4,21].At g x = 0 we find that the Hopf bifurcation ends at the lower saddle-node type bifurcation surface at g x = 0 in a codimension-2 Takens-Bogdanov line.On this line the Jacobian has a double zero eigenvalue.In addition to the saddle-node type bifurcation and the Hopf bifurcation a homoclinic bifurcation emerges from the Takens-Bogdanov bifurcation line.This bifurcation is in general difficult to detect and can ecologically be related to sudden population bursts.In the model class under consideration a Takens-Bogdanov bifurcation can only be observed in systems with non monotonic functional response G(X).This property is necessary to enable negative values of g x and therefore it is also necessary to cross the Takens-Bogdanov bifurcation at g x = 0.
Another way to achieve a destabilization of the steady state due to a Hopf bifurcation is to decrease the competition parameter m x .This effect is rather counterintuitive since decreasing the competition means ecologically to improve the food conditions for the predators.Such behavior of the model can be strongly related to the paradox of enrichment [47].

Disease in the predator population
We now investigate the impact of a disease in this two trophic food chain as described by Eq. ( 3).In contrast to the normalized predator-prey model the normalized eco-epidemic model has not 3 but 13 parameters that may all more or less influence the stability of the steady state.Thoroughly analyzing and discussing these parameter variations is beyond the scope of this work.Instead we focus only on bifurcations that give evidence for more complex dynamics.Such complex behavior like quasiperiodic or chaotic dynamics occur usually in the neighborhood of global bifurcations or bifurcations of higher codimension, which can be found for a lot of different parameter sets in the model class under consideration.It is important to note that though the whole bifurcation analysis presented is based on local stability properties of the steady state, we are able to detect easily higher codimension bifurcations since they correspond to intersections of different bifurcation surfaces like the Takens-Bogdanov line discussed in the previous subsection.
For our analysis we chose g x , l s and l i as the most important bifurcation parameters.For most of the common incidence functions (cf.Table 1) l i ≤ 1 due to saturation effects with respect to the number of infected.The sensitivity to prey g x is for the same reason also confined to this range g x ≤ 1.
For the the scaling parameters we assume now that 90 percent of the losses of prey are caused by predation, which means a relative competition m x = 0.1.We assume further that 95 percent of the predation is caused by healthy predators, i.e. susceptibles plus recovered, ( fα = 1 − f α = 0.95) and that half of the healthy predators are susceptible (b = 0.5).The gain of susceptible predators results for 95 percent from biomass conversion (e s = 0.95) and only 5 percent from the recovering (ẽ s = 1 − e s = 0.05).The natural mortality of the healthy predators is assumed to be relatively low compared to losses due to infection (m y = 0.1).
Using the parameter settings mentioned above we compute the stability of the steady state as shown in Fig. 2. We find as in the model without disease a surface of saddle-node type bifurcations (transparent blue) and a surface of Hopf bifurcations (red).But now the Hopf bifurcation surface possesses a rather complicated shape.This shape corresponds to a Whitney umbrella, a bifurcation situation which is rarely found in applications.Other examples of a Whitney umbrella are presented in [23,52].In this bifurcation scenario the Hopf bifurcation surface is twisted around a codimension-3 1:1 resonant double-Hopf point characterized by two identical pairs of complex conjugate eigenvalues.As Fig. 2 shows, a line of codimension-2 double-Hopf bifurcations emerges from this point.At this line where the Hopf bifurcation surface intersects itself two pairs of purely imaginary complex conjugate eigenvalues can be found.
Additionally we observe two intersection lines of the Hopf bifurcation surface with the saddle-node type bifurcation surface.One is again a Takens-Bogdanov bifurcation and the other one is a Gavrilov-Guckenheimer bifurcation line that emerges from a so-called triple point bifurcation on the Takens-Bogdanov bifurcation line [38].On this Gavrilov-Guckenheimer line the Jacobian has a zero eigenvalue in addition to a pair of purely imaginary complex conjugate eigenvalues.In contrast to the Takens-Bogdanov bifurcation the Hopf bifurcation surface does not end on the Gavrilov-Guckenheimer bifurcation.The existence of the Gavrilov-Guckenheimer bifurcation indicates that quasiperiodic and chaotic dynamics are likely to occur in the neighborhood of this bifurcation.The double-Hopf bifurcation line instead is a clear evidence for the existence of chaotic parameter regions [38].Therefore we can conclude that the consideration of the vulnerability of the predator population to a disease can lead in general to complex dynamics in eco-epidemiological systems.
Unfortunately we have no information about the size of the chaotic parameter region since our analysis is based on mathematical theorems.In the following section we investigate a specific model that allows to translate the generalized parameters to specific system parameters and vice versa.This specific example allows not only to explicitly compute the chaotic parameter regions but also gives insights into the route to chaos.

Chaos in a specific eco-epidemiological system
To demonstrate the theoretical implications of the existence of a Whitney umbrella bifurcation situation in a specific model, we need to choose specific mathematical functions for the generalized processes.This means that we construct an example of a specific model that compares to the generalized parameter set of the bifurcation diagram shown in Fig. 2.
Since the generalized parameters represent the derivatives of the generalized processes we need to choose functions according to the parameter range of the higher codimension bifurcations.As stated above in ecology and epidemiology a large pool of proposed functional responses G(X) and incidence functions λ(Y S , Y I ) exists.Since the double-Hopf bifurcation line appears for g x lower than 1, which means a increase of g(x) slower than linear in x, a functional response with saturation is necessary to obtain the double-Hopf bifurcation line.
Instead of defining a specific model and normalizing it we construct an already normalized specific model for the sake of simplicity.We choose a Holling type III function g(x) = ax 2 /(1 + bx 2 ) as the functional response.Due to the normalization we need g(1) = 1.Therefore we define a := (1 + b).The relation between b and g x is then g x = 2/(1 + b).In a similar way we allow values of l s and l i lower than 1 as well.We use the asymptotic incidence function l(y s , y i ) = cy s y i /(1 + dy s + ey i ) [17,46].In order to satisfy l(1, 1) = 1 we define c := (1 + d + e).We find the relations l s = 1 + e/(1 + d + e) and Now we have a specific model that allows a translation of the parameter set of the generalized model into the parameters of the specific model.This enables us to analyze the dynamics of the system.As discussed in the previous section the double-Hopf bifurcation indicates the emergence of chaotic parameter regions.In order to find these parameter regions we compute the Lyapunov exponents of the specific system for a grid of points in the generalized parameter space close to the double-Hopf bifurcation.
The result is shown in Fig. 3.The two solid lines are Hopf bifurcation lines that intersect in a double-Hopf bifurcation.Within the white area the steady state (1,1) is stable.In the light grey area the system exhibits periodic long-term dynamics.
In addition to the two Hopf bifurcation lines we find numerically other bifurcation lines (two dashed, one dotted) using pathfollowing methods implemented in MATCONT [16].The dashed lines are Neimark-Sacker bifurcations where a limit cycle becomes unstable and a stable quasiperiodic motion on a torus emerges.This behavior can be found in the dark grey parameter regions.In the black regions the largest Lyapunov exponent is positive and hence, the dynamics is chaotic.The dotted line is a period doubling bifurcation.For small values of l s we find first the transition to quasiperiodic motion on a torus with a subsequent transition to chaos.For larger values of l s the periodic solution undergoes first a period doubling before the torus or Neimark-Sacker bifurca- Beyond the Hopf bifurcation lines (1.) the system exhibits stable periodic dynamics (bright grey).Both Hopf bifurcation lines intersect in a codimension-2 double-Hopf bifurcation.This double-Hopf bifurcation is the starting point of a Neimark-Sacker bifurcation line (dashed).At the Neimark-Sacker bifurcation line a stable Torus emerges (2.) and the system exhibits quasiperiodic dynamics (dark grey).The dotted line is a period doubling bifurcation line.Beyond this line where the system oscillates on a period 2 orbit (3.) we find another Neimark-Sacker bifurcation (dotted-dashed line).In the black region we find chaotic dynamics (4.,5.).
tion occurs.While the transition to chaos involves always a transition from quasiperiodicity, the chaotic attractor looks different for small and large values of l s since the Neimark-Sacker bifurcation and the period doubling swap places.Both routes to chaos are illustrated in Fig. 4.
The population dynamical system alone (Eq.( 1)) as well as the epidemiological system alone (Eq.(2)) do not exhibit complex dynamics.Only if both are coupled to form an eco-epidemiological system with an infected predator, the dynamics can be quasiperiodic or chaotic.Our generalized analysis shows that chaos is generic in this class of models.), a torus at g x = 0.6025, l s = 0.456 (2.), a limit cycle with doubled period at g x = 0.613, l s = 0.5 (3.) and two chaotic attractors at g x = 0.6, l s = 0.47 (4.) and g x = 0.608, l s = 0.51 (5.).

Discussion
We have studied a generalized eco-epidemic model which couples the behavior of a predator-prey system to the dynamics of a disease which can infect the predator.The advantage of investigating generalized models lies in the fact that the exact mathematical form of the interaction processes like predatorprey or infection interactions does not have to be specified.This allows for rather general conclusions about the stability of the positive steady state which will be reached in the long-term limit.Moreover, this generalized approach can give insight into the global dynamics of the system though only a local stability analysis is performed.Due to the usage of generalized models our results apply to certain classes of models.
The eco-epidemic model is based on two often used generalized models which exhibit only stationary points or periodic behavior when studied separately.
We have shown that the coupling of an ecological and an epidemiological model can lead to classes of systems exhibiting complex dynamics like quasiperiodic and chaotic behavior.Our result is based on the detection of higher codi-mension bifurcations like double-Hopf bifurcations, triple points and Gavrilov-Guckenheimer bifurcations.In the neighborhood of such bifurcations there exist parameter regions where quasiperiodic and chaotic behavior can be found.This mathematical finding based on bifurcation theory is illustrated with a specific model where the interaction functions are specified in order to be able to use numerical methods like path-following of bifurcations and the computation of Lyapunov exponents.As a result we can demonstrate that the chaotic parameter regions are not small and therefore not negligible, but rather large and hence, important for the dynamics of the system.
Finally we note that we were not able to find complex dynamics in the ecoepidemic model using a SIS instead of a SIRS model.Therefore it seems to be essential that we introduce the class of recovered predators for the occurrence of complex dynamics.
To our knowledge this is the first example where chaotic behavior has been found in an eco-epidemic model system with a disease in the predator population.Based on the generalized approach we can state that the emergence of chaos is generic for certain classes of eco-epidemic models and thus likely to be found.

g x sensitivity toFig. 1 .
Fig. 1.Bifurcation diagram of a generalized predator-prey model.A surface of Hopf bifurcations (dark) and two surfaces of saddle-node type bifurcations (transparent bright) are shown.The bifurcation parameters are the prey sensitivity g x , the timescale of the predator a y and the competition m x (intraspecific competition of the prey).

Fig. 2 .
Fig. 2. Bifurcation diagram of the generalized eco-epidemic model.A surface of Hopf bifurcations (dark) and a surface of saddle-node type bifurcations (transparent bright) are shown.The intersection lines are a Takens-Bogdanov bifurcation line (TB), a Gavrilov-Guckenheimer bifurcation line (GG) and a double-Hopf bifurcation line (DH).The double-Hopf bifurcation line ends in a 1:1 resonant double-Hopf bifurcation point (1:1 DH) and the Gavrilov-Guckenheimer bifurcation line ends in a triple point bifurcation at the Takens-Bogdanov line.The bifurcation parameters are the generalized parameters g x , l s and l i which are strongly related to the functional form of the underlying processes, namely the per capita functional responce g(x) and the incidence function l(y s , y i ) respectively.The fixed parameters are scale parameters a x = 1, a s = 0.4, a i = 0.9, a r = 0.25, f α = 0.05, b = 0.5, e s = 0.98, m x = 0.1 and m y = 0.1.

Fig. 3 .
Fig. 3. Dynamics of a specific model close to the double-Hopf bifurcation.The steady state is stable in the white area between the two Hopf bifurcation lines (solid lines).Beyond the Hopf bifurcation lines (1.) the system exhibits stable periodic dynamics (bright grey).Both Hopf bifurcation lines intersect in a codimension-2 double-Hopf bifurcation.This double-Hopf bifurcation is the starting point of a Neimark-Sacker bifurcation line (dashed).At the Neimark-Sacker bifurcation line a stable Torus emerges (2.) and the system exhibits quasiperiodic dynamics (dark grey).The dotted line is a period doubling bifurcation line.Beyond this line where the system oscillates on a period 2 orbit (3.) we find another Neimark-Sacker bifurcation (dotted-dashed line).In the black region we find chaotic dynamics (4.,5.).