A Mathematical Model Supporting a Hyperpredation Effect in the Apparent Competition Between Invasive Eastern Cottontail and Native European Hare

In this work a mathematical model is built in order to validate on theoretical grounds field study results on a three-species system made of two prey, of which one is native and another one invasive, together with a native predator. Specifically, our results mathematically describe the negative effect on the native European hare after the introduction of the invasive Eastern cottontail, mediated by an increased predation rate by foxes. Two nonexclusive assumptions can be made: an increase in cottontail abundance would lead to a larger fox population, magnifying their predatory impact (“hyperpredation”) on hares; alternatively, cottontails attract foxes in patches where they live, which are also important resting sites for hares and consequently the increased presence of foxes results in a higher predation rates on hares. The model results support hyperpredation of increasing fox populations on native hares.


Introduction
Biological invasions-i.e. the human-mediated introduction of species outside their native range-represent one of the main drives of global change. Introduced species change ecosystems composition and species interaction, threatening native biodiversity and representing a major source of extinction (Clavero and García-Berthou 2005;Kumschick et al. 2015). Introduced species may interact negatively with native ones, through numerous, such as competition, predation, hybridization, the transmission of disease, see IUCN 2020 for a complete list of impact mechanisms (IUCN EICAT Categories and Criteria 2020). Competition is the interaction between species of the same trophic level, which leads to negative consequences for one of the two. It can be regarded as the direct or indirect interaction of organism that leads to a change in fitness when they share the same resources, such as food or nesting sites (Holomuzki et al. 2010). The reduction in the fitness of individuals translates at the population level into a change in demography which could determine the decline or the extinction of one of the two species. There are three major mechanisms of competition. Interference competition and exploitation competition are categorized as real competition. Interference competition involves direct interactions between individuals. For example, in Ferretti et al. (2011) it is reported the negative effects of behavioural interference by fallow deer (Dama dama) on the foraging behaviour of roe deer (Capreolus capreolus). Exploitation competition implies indirect negative interactions between two species arising from the use of a common resource (Schoener 1983). A classic example is the competition between the introduced Eastern grey squirrel (Sciurus carolinesis) and the native Eurasian red squirrel (Sciurus vulgaris). The replacement is mainly due to exploitation competition for food resources, with the introduced species more efficient in their use (Wauters et al. 2002a, b). The third form of competition is more complex and concerns species that do not interact directly and do not exploit the same resources, but still influence each other through shared enemies, such as predators, parasites, or pathogens (Holt and Bonsall 2017). In this case, the competition is mediated by the action of a third species of another trophic level. This species might be a pathogen: in Great Britain, grey squirrels carry a squirrel pox virus which is lethal for the red squirrel (Tompkins et al. 2003;Romeo et al. 2019). The virus carried by the introduced species has a detrimental effect on the native one, increasing the replacement speed between the two species. However, the common enemy might also be a predator: in New Zealand, introduced rabbits (Oryctolagus cuniculus) created large populations of mammal predators, which also prey upon native lizards (Oligosoma spp. (Norbury 2001).
The Eastern cottontail (Sylvilagus floridanus) is a lagomorph native to the American continent that was introduced in Italy in the 1960s for hunting purposes (Bertolino et al. 2011a). The species is now widespread in Central and Northern Italy (Bertolino et al. 2011a;Loy et al. 2019) where it competes with the native European hare (Lepus europaeus). The two species select different macro-and microhabitats, both for resting and during feeding activity (Bertolino et al. 2011a(Bertolino et al. , 2013. Cottontails carry several viruses and parasites, which can potentially affect hares (Bertolino et al. 2010;Tizzani et al. 2014); however, competitive interactions mediated by parasites have not yet been highlighted. Recently, in Cerri et al. (2017) an apparent competition between Eastern cottontail and European hare mediated by the fox (Vulpes vulpes) as a common predator is shown. Introduced cottontails affect prey-predator dynamics of native hares and foxes: the abundance of foxes is positively associated with the one of hares, when cottontails are scarce, suggesting the main influence of external factors, such as habitat quality. However, this association becomes increasingly negative as cottontails increase in their abundance. When cottontails are abundant, increases in the abundance of foxes are negatively correlated with the abundance of European hares (see Fig. 2 in Cerri et al. (2017)). This pattern is suggestive of an indirect competition between cottontails and hares mediated by fox predation.
Starting from the work of Cerri et al. (2017), we introduce here a mathematical model to simulate this three-species system and investigate it, to possibly validate on theoretical grounds the results of the field study (Cerri et al. 2017).
The paper is organized as follows. After the model formulation in the next section, its equilibria are assessed in Sect. 3. The three-species system behaviour is analysed in detail in Sect. 4, numerical simulations are reported in Sect. 5, and a biological discussion sets these results in the ecological perspective, in Sect. 6.

Model Formulation
The dynamical system that we propose can be described by the following equations, where all the parameters are nonnegative and where the parameter r generally interpreted as the traditional growth rate of the logistic equation, represents instead the foxes reproduction. This important difference should be remarked throughout the paper.
The reason for separating the reproduction r and the mortality m rates for the foxes, and more in general for the various model populations, and not using the standard logistic equation, lies in the fact that we are going to use field data, for which, in principle, the condition of the positivity of the net reproduction rate may not be satisfied, and furthermore, that we would like to contemplate the possibility of disappearing populations.
We also use the Holling type (HT) I response function for hunting because the density of the foxes is very low, while hares and cottontails do not attain very large numbers. Indeed, although not much is known about the foxes density in Italy, a reasonable range per square kilometre is usually taken as [1.0, 2.5], (Boitani and Prigioni 2003). The maximal observed cottontail density is 110 per km 2 (Bertolino et al. 2011a), while for hares it ranges in the interval [26,40], (Pandini et al. 1998). Therefore, a feeding saturation phenomenon for foxes does not arise, and a simple bilinear term is sufficient to describe the interactions, as it represents a linear approximation in the meaningful range of the HTII response function describing satiation.
The first equation in (1) describes the dynamics of the red fox Vulpes vulpes V that grows logistically with reproduction rate r , death rate m, intraspecific competition rate c V V and gets a benefit from the capture of its prey scaled by a conversion coefficient e, related to the digestibility of the prey, which, for foxes, lies in the range [0.8958, 0.9190], when predating on rabbits (Maurya et al. 2012), while hunting hares, it is 0.910 (Ruehe et al. 2008). As these figures are similar, we take the conversion factor e to be the same for the two kinds of prey.
The second equation describes the dynamics of the Eastern cottontail Silvilagus floridanus, S, that also grows logistically with reproduction rate s, mortality rate n, has the intraspecific competition rate c SS and is hunted by the foxes at rate a.
The last equation describes the dynamics of the European hare Lepus europaeus, L. The logistic growth occurs with reproduction rate u, mortality rate p, intraspecific competition rate c L L . Hares are captured by foxes at rate b. In Table 1 we summarize the meaning of each model parameter.
Note that no direct competition between the invasive cottontails and the native hares is here considered, as on biological grounds the two species occupy different ecological niches and thus do not directly interfere with each other (Bertolino et al. 2011a(Bertolino et al. , b, 2013. Rather, it is hypothesized that their interaction occurs through indirect competition mediated by the red fox.
The model (1) contains the predators' interference on the two types of prey, but the latter occupy different ecological niches (Bertolino et al. 2011a(Bertolino et al. , b, 2013, for which, in the absence of the foxes, each one individually or both together settle to their respective carrying capacities, that are explicitly given by the nonvanishing population values contained in the equilibria E 1 , E 2 , E 4 and E 6 below. When foxes are present and one prey is absent, the population levels attained are clearly modified by their mutual interactions, see equilibria E 3 , E 5 . For the later study of the equilibria stability, we will need the Jacobian of the system (1):

Boundedness
In this subsection we show that the system trajectories are confined in a compact set in the first quadrant of the phase space. This result is relevant from the ecological point of view, because it says that in the presence of finite resources, no population can grow without limit. The result is obtained by considering all the living individuals in the model. Let us define the total environmental population A = V + S + L. Because each population is nonnegative, if A is bounded, then also V , S and L are bounded. The details for which boundedness of A holds are reported in "Appendix".

Equilibria
System (1) has the following eight possible equilibria The origin, E 0 = (0, 0, 0), corresponds to extinction of the whole threespecies system. The three points in which only one population thrives: The equilibria with only one vanishing population: and coexistence E 7 = (V 7 , S 7 , L 7 ) whose population levels can explicitly be evaluated: While the origin is always feasible, it turns out to be conditionally stable, in view of the fact that the Jacobian reduces to a diagonal matrix from which the eigenvalues are immediate and provide the following stability conditions: (2) Note that these conditions imply indeed that mortalities exceed the reproduction rates, entailing the extinction of each and every species. For E i , i = 1, 2, 4, the feasibility conditions are, respectively, r ≥ m; s ≥ n; u ≥ p.
For E 1 stability hinges on the following inequalities, m < r , which is implied by feasibility and is therefore redundant, and: For E 2 we find instead the stability conditions and n < s, which follows from feasibility. For E 4 , which is feasible if the hares population is nonnegative, i.e. the last condition in (3) that ensures also the negativity of the last eigenvalue, the stability conditions come from the other two explicit eigenvalues and give: For the equilibria with two nonvanishing populations, the results are as follows. E 3 is feasible for Here one eigenvalue factorizes, namely J 33 (E 3 ) = u − p − V 3 . The remaining submatrix has negative trace and positive determinant; hence, the Routh-Hurwitz conditions are satisfied and thus, its eigenvalues are negative. Stability reduces to requiring At E 5 two populations need to be nonnegative, giving for feasibility As for E 3 , here one eigenvalue can be factorized, J 22 (E 5 ) = s−n−aV 5 . The remaining minor has negative trace and positive determinant, and therefore, its eigenvalues are negative. The only stability condition is: For E 6 feasibility follows by satisfying and all the Jacobian eigenvalues are explicitly known, r −m+eaS 6 +ebL 6 , −sc SS < 0, −uc L L < 0. Thus, just the first one ensures stability, namely: Feasibility for E 7 entails the following inequalities: This equilibrium, when feasible, turns out to be unconditionally stable. Indeed, using the Routh-Hurwitz criterion, for the trace we find: and for the determinant: Further, we need to assess the sum M of the principal minors of order 2: The final Routh-Hurwitz condition for stability is then: This condition is always satisfied, as it reduces to

Equilibria Global Stability
The conditions summarized in Table 2 suggest that there are pairs of equilibria that cannot occur simultaneously. This will be better investigated in the next subsection, but since bi-or multistability among some subsets of equilibria is prevented by these conditions, it is worth to investigate their possible global stability. To this end, we construct suitable Lyapunov function candidates L i (V , S, L) for each equilibrium, such that nonnegativity, in the first orthant, as well as It turns out that all equilibria are globally asymptotically stable, whenever they are locally asymptotically stable. In fact it turns out that the matrix A i is independent of the equilibrium that is considered, i = 1, . . . , 7; it has namely the following structure: The mathematical details are deferred to "Appendix".

Bifurcations
In view of the results of the previous subsection, Hopf bifurcations are forbidden, a fact that prevents the onset of persistent oscillations in the model solution trajectories. Indeed, at every equilibrium, either the eigenvalues are all explicitly known and real, or, when the Routh-Hurwitz conditions are used, in case of equilibria E 3 and E 5 , the trace can never vanish, and this prevents the possibility of having pure imaginary eigenvalues. In case of E 7 instead, it is the condition (14) that cannot vanish. In saying so, we of course exclude, as biologically unrealistic, the very particular cases in which some or all the reproduction rates of the three species vanish. Therefore, only transcritical bifurcations could then relate the various system's equilibria to each other. We will now rigorously explore this issue. The mathematical tool is represented by Sotomayor's theorem, (Perko 2001).
, we first evaluate its partial derivatives with respect to the model parameters: We then need the Jacobians of the above vector-valued functions: Then, in order to evaluate D 2 F we need: Further, to evaluate D 3 F we need the partial derivatives of order three, but it can easily be assessed that they all vanish. As a consequence, system (1) is therefore not satisfying the necessary condition for a pitchfork bifurcation: Therefore, all pitchfork bifurcations are ruled out and so possibly only saddle-node or transcritical bifurcations can arise. The mathematical details for this analysis are reported in the "Appendix", and as we will see, only the latter are indeed found. A picture of the structure of all the analytically found bifurcations is reported in Fig. 1.

Biological Interpretation
We now provide some insights on the situation, based on all the previous findings, in terms of some relevant biological quantities. In particular, we will focus on the net reproduction rates of each species and some additional quantities, namely: In addition, recall that bistability is forbidden. Thus, in the following discussion, as always the equilibria with some of the vanishing populations appear, coexistence is impossible for all the listed situations.
-When all the net reproduction rates are negative r < m, s < n, u < p only E 0 is feasible and stable; therefore, as intuition suggests, all the species are driven to extinction. -For the foxes reproducing effectively only E 1 is feasible and stable: This makes sense, as foxes are supposed to have other feeding resources.
-When r < m, n < s, u < p E 2 is feasible, but stable only if c SS < h 2 , while in the opposite case c SS > h 2 , E 3 is feasible and stable. In this situation, the cottontails thrive alone if their intraspecific competition is low enough, quite intuitively; otherwise, they would support also the foxes and the two species survive together.
-Instead, if m < r , n < s, u < p we find only the equilibria E 1 and E 3 . The former is feasible and stable only if c V V < k 1 , the latter for feasibility needs c V V > k 1 . Thus, here either the foxes survive alone, if their intraspecific competition is low enough; otherwise, they coexist with the cottontails. -For r < m, s < n, p < u the possible equilibria are E 4 and E 5 . E 4 is unconditionally feasible, but it is stable only if c L L < g 2 . Conversely, for c L L > g 2 , E 5 becomes feasible and in such case it is stable. Thus, here either the hares survive alone, provided their intraspecific competition falls below a threshold; otherwise, they thrive together with the generalist predators, the foxes. -When m < r , s < n, p < u E 1 is feasible but stable only if c V V < k 2 , while for the opposite condition, E 5 is feasible and always stable. Here for a small enough intraspecific competition, the foxes thrive alone, wiping out both hares and cottontails, but in the opposite case the hares also survive. -Whenever r < m, n < s, p < u foxes and cottontails coexist at E 3 , or the foxes thrive with the hares at E 5 or else foxes disappear and only cottontails and hares share the environment. More precisely, E 3 is feasible if c SS < h 2 and stable if u < p + bV 3 , E 5 is feasible if c L L < g 2 and stable for s < n + aV 5 , E 6 is unconditionally feasible but stable only if r + eaS 6 + ebL 6 < m.
-Finally, for m < r , n < s, p < u only one population survives unconditionally in the environment, the foxes; here all possible alternative equilibria are feasible, and only their stability determines the three-species interaction outcome. More precisely, the foxes wipe out the other two populations, equilibrium E 1 , if c V V < k 1 and c V V < k 2 , they thrive with the cottontails at E 3 if c SS > h 2 , c V V < k 1 and u < p + bV 3 , these conditions, respectively, arising for feasibility and stability, or finally, they survive together with the hares, equilibrium E 5 which is feasible for c V V < k 2 and c L L < g 2 , while stable in case s < n + aV 5 . To better illustrate how the species at steady state are related to each other, we provide Fig. 1, where the graph contains as nodes the equilibria, whose nonvanishing populations are denoted by subscripts, and the arcs their connecting transcritical bifurcations.

Numerical Simulations
Although the demography of three-species system has been completely characterized analytically, it is worth to see also some simulations with realistic parameter values, in order to compare their results with the field findings of the biologists.

Parameter Values
Some of the parameter values are taken from as well as from (Amori et al. 2008;Barbara et al. 2018;La Morgia and Venturino 2017).
For mortalities, the average lifetimes are taken from (Amori et al. 2008), and give for Lepus europaeus a range of 5-6 years, for Sylvilagus floridanus it is about 15 months, with individuals living up to 5 years, and for Vulpes vulpes it seldom exceeds 3-4 years. We therefore take the following averages m = 1 3.5 = 0.28571, n = 1 1.25 = 0.80000, p = 1 5.5 = 0.18182.
For the conversion coefficient e, i.e. the benefit that foxes obtain hunting their prey, as mentioned in Introduction, we set e = 0.91.
To assess the reproduction rate r P for a generic population P, we consider the following equation:Ṗ = r P P whose solution is: P(t) = e r P t P(0). Now, each female fox gives to life an average of 4 offsprings so that the foxes population will triplicate in 1 year: P(1) = 3P(0). The foxes population growth is thus obtained solving 3P(0) = e r P(0) for r and obtaining r = log 3. Data on the densities of fox populations in Italy are scanty. They probably range between 1.0-2.5 foxes per km 2 (Boitani and Prigioni 2003). Here we use a conservative value of 1.0 fox per km 2 . Using (16) and solving for the intraspecific competition term from the equilibrium E 1 , we find For the cottontails, the biologist data indicate that the cottontails population in Northern Italy after the reproductive period increases 4.5 times; this of course already discounts mortality, so that we can take P(1) = 4.5P(0), giving s = log 4.5. The maximum density recorded in lowland in Piedmont, the same region where the study of Cerri et al. (2017) was performed, is 110 cottontails per km 2 (Bertolino et al. 2011a). Since this value was only from an area, we use instead 100 cottontails per km 2 . We then find the intraspecific competition coefficient from the equilibrium E 2 and (16), c SS = 1 100 log 4.5 − 1 1.25 = 0.0070408.
The hares reproduction is highly variable; each female can give birth to 1-4 (with a mean of 2.6) juveniles and reproduces 2-5 times a year. We considered a mean of 2.6 juveniles with 3 reproductions and an annual output of 8 young ones. Hence, the population quintuplicates in 1 year, P(1) = 5P(0); we thus set the reproduction as u = log 5. The carrying capacity for the European hare is 26-40 hares per km 2 in territories with a good suitability (Pandini et al. 1998), i.e. mixed crops with wheat, meadows and other cultures. Taking 30 hares per km 2 and combining E 4 with (16), we thus find the intraspecific competition rate c L L = 1 30 log 5 − 1 5.5 = 0.047587.

The Bifurcations Chains
We provide now a few figures that relate some of the transcritical bifurcations to each other, where the transitions between one equilibrium and the next one depend on one specific model parameter. The parameter reference values are those of Table 3. For each case instead we provide the remaining parameters.
Taking as bifurcation parameter a, the hunting rate on cottontails, for its increasing values we find three bifurcations, located at In turn, the three-species interaction moves from coexistence, to the point where hares vanish, to coexistence again, and finally to the foxes-hares subsystem, see Fig. 2.
The transition E 5 − E 1 Here the bifurcation parameter is b, the hunting rate on hares, in the range [16,22]. For its increasing values, the three-species system moves from coexistence of both foxes and hares to the situation in which only foxes thrive, in view of the possibility of using other food supplies, Fig. 3. The specific parameter values are: The transition E 3 − E 1 Taking again a as an increasing bifurcation parameter, the three-species system moves from foxes-cottontails coexistence, to the point at which only foxes thrive, Figure Online) For small values of a, we find coexistence, equilibrium E 7 , then at a 1 3,7 = 1.0078 the hares vanish giving equilibrium E 3 , subsequently in the interval [a 2 3,7 = 2.0762, a 5,7 = 4.8392] hares reappear, and coexistence is re-established. Finally, for values of a exceeding a 5,7 = 4.8392, cottontails are wiped out and the system attains equilibrium E 5 (Color Figure  Online

Fig. 2 (Color
The transition E 7 − E 5 − E 7 − E 3 Figure Online) For b small, coexistence is obtained, equilibrium E 7 . At b 1 5,7 = 1.6581 cottontails disappear, leaving hares and foxes in the environment, then at b 2 5,7 = 3.6921 they reappear, so that coexistence of the three species is re-established, finally at b 7,3 = 5.7296 the hares are wiped out, and cottontails thrive together with foxes (Color Figure Online for which the three-species system moves from coexistence to the hares-foxes subsystem, to coexistence again, and finally to the cottontails-foxes subsystem. The remaining parameters here are:

Fig. 5 (Color
The transition E 5 − E 7 − E 3 The bifurcation parameter b, the predation rate on hares, in the interval [3, 7] generates the bifurcations for which the system undergoes a transition from E 5 , the cottontails-free equilibrium, to E 7 , coexistence, and finally to E 3 , the hares-free point, Fig. 6. The other parameters are: Note that this trend is biologically interesting, because it reflects a possible behaviour observed in the field studies. Indeed, in the areas under observation, the cottontails are an invasive species that modified the natural predator-prey dynamics. The results of the field study (Cerri et al. 2017) indicate that a larger cottontails population drives hares to extinction. As it can be seen from Fig. 6, increasing the predation rate on hares, the cottontails population increases and finally replaces the hares, because the latter are subject to an increasing hunting rate by the foxes, and are therefore driven to extinction.

Discussion
The question raised in Introduction has therefore an answer. the invasion mechanism becomes apparent. The cottontail population increases steadily, and following it, also the foxes attain higher population values. In turn, their damage to the hares increases, with the latter dropping and finally becoming extinct, see Fig. 7. In particular in the figure one can observe that while the cottontail population invades, growing from low values at the early times to a larger population around time t = 0.4, the foxes population also raises up. Thus, their pressure on their prey increases, but the hares suffer more, their numbers decrease and eventually they are wiped out. Note that this happens even though the predation rate on cottontails exceeds the one on hares, b = 3.5 > a = 1.5.
To further and strengthen our analysis, we try to relate our findings to the ones coming from field observations in Cerri et al. (2017). In particular, we found a relationship between hares and foxes similar to the one described in the last three figures of Cerri et al. (2017), by running a set of numerical simulations with randomly generated 51 Page 18 of 28 E. Caudera et al. Figure Online) Logarithm of hares population versus foxes. Left to right and top to bottom, each frame is related to higher values of cottontails densities. Specifically, top: 1 < S < 3, 3 < S < 5, 5 < S < 7; bottom: 7 < S < 9, 9 < S < 11, 11 < S < 13 (Color Figure Online) parameter values, so as to mimic the outcome of the field data gathered over an 8 year timespan.

Fig. 8 (Color
We thus ran 500 simulations of the model (1), taking uniform random values for all the parameters, located in intervals of semiwidth 9% above and below the reference values of Table 3 supplemented by the values indicated here: The steady states were then evaluated and recorded for each simulation, together with the respective parameter values that generated them. Finally, regression lines were evaluated, and the results plotted again for the logarithm of the hares population versus the foxes. Two selected such results are reported in Figs. 8, 9, corresponding to two different random parameter choices. The various frames, left to right, report the findings for increasing values of the cottontails density. The plots represent the logarithm of the hares populations versus the foxes. There is a remarkable resemblance with the analogous plot obtained by the field data and reported in Cerri et al. (2017).
The transition E 7 − E 3 − E 7 − E 5 shows a situation where the hunting rate of foxes on cottontails is not limited and continues to increase over time. In this case, the system moves from an initial situation where predation of foxes on the introduced species is limited and the three species coexist to a final outcome where cottontails are wiped out and only the two native species remain. Such a situation would be positive for biodiversity because a native predator would thrive to extinction the introduced lagomorph. However, the red fox is a generalist predator, (Soe et al. 2017), and although the cottontails became an important prey of foxes in Italy (Balestrieri et al. 2006), the availability of many other food resources (e.g. rodents, insects, fruits) probably prevents the predator from increasing the predation rate beyond a certain limit. Left to right and top to bottom, each frame is related to higher values of cottontails densities. Specifically, top: 1 < S < 3, 3 < S < 5, 5 < S < 7; bottom: 7 < S < 9, 9 < S < 11, 11 < S < 13 (Color Figure Online) Conversely with the transition E 7 − E 5 − E 7 − E 3 , the hare declines rapidly starting from a very high and unrealistic density. The cottontail goes extinct and then recovers when the hare is at low density. This would indicate a strong predation pressure on the hare and a competitive effect of hares on cottontails, which can recolonize the area only when the hare is at low density. It is known that the fox cannot limit hares so drastically and that hares do not exclude cottontails (Bertolino et al. 2011a).
The transitions E 5 − E 1 and E 3 − E 1 indicate that the fox predation pressure on cottontails and hares is so high as to bring the two prey species to extinction. This situation is not supported by earlier field research. In areas where only the two native species are present, the hare is one of the food sources for foxes. However, predation is not considered to be a limiting factor for the lagomorph (Boitani and Prigioni 2003;Cerri et al. 2017). Indeed, in places where the hare has disappeared and the cottontail is the only widespread lagomorph, the American species becomes a main prey of foxes, but populations survive due to their high reproduction rates (Balestrieri et al. 2006).
The transition E 5 −E 7 −E 3 better reflects the results of previous research (Cerri et al. 2017). The system undergoes a transition from an alien species-free system to a one invaded by the introduced cottontail; subsequently, it evolves to a coexistence regime between the three species and finally the hare becomes extinct. For the parameter values (22), the simulations of this model show an agreement with the invasion process recorded on the field (Cerri et al. 2017). The simulation in Fig. 7 describes the invasion of an area by cottontail with a rapid increase in its population and a benefit for foxes that exploit this new prey and in turn increase their density. The predation pressure on both the hares and the cottontails consequently increases, but the hares suffer more and eventually go extinct. The cottontail better supports the higher foxes predation rate, due to its higher reproductive output (Chapman et al. 1980).

Conclusion
In this work we were able to simulate a three-species system composed of two prey species, one native and one introduced, and their relations with a native predator. Our results mathematically describe the negative effect on the native European hare after the introduction of the invasive Eastern cottontail, mediated by an increased predation rate by foxes (Cerri et al. 2017). In Cerri et al. (2017) two nonexclusive hypotheses were discussed to explain the pattern observed when cottontail populations increase. In the first hypothesis, an increase in cottontail abundance would lead to a numerical response of foxes, magnifying their predatory impact on hares (a sort of "hyperpredation", a particular case of apparent competition). The second hypothesis assumes that cottontails attract foxes in patches of permanent cover where they live. Since these habitats are also important resting sites for adult and young hares (Boitani and Prigioni 2003), an increased presence of foxes is likely to result in a higher predation risk for hares. Our simulation describes an increase in foxes' densities after the invasion of cottontail, supporting the first hypothesis from Cerri et al. (2017) based on a hyperpredation of increasing fox populations on native hares.
Thus, the upper bound is found for which Upon integration of this differential inequality, the final desired bound is obtained:

Global stability details
Equilibrium E 0 In this case we set where c, g and f are arbitrary positive coefficients. Differentiating L 0 with respect to time, we obtain, recalling (15), where the last inequality follows on using the stability conditions (2). Now, choosing f = g = e and c = 1, A becomes diagonal, with negative eigenvalues; hence, it is negative definite as required.
Equilibrium E 1 Here we set Upon differentiation, we find having used the stability conditions (4) in the last step. Choosing f = g = e and c = 1, A is seen to have negative eigenvalues, and in turn is negative definite. Global stability of E 1 therefore follows.
Equilibrium E 2 The candidate function is which follows using (5). Negative definiteness is then obtained by the choice f = g = e and c = 1.
Equilibrium E 3 We now set Differentiating with the help of (8), we obtain Choosing f = g = e and c = 1, A turns to be negative definite. Equilibrium E 4 We now choose which, together with (6), gives Global stability follows by setting f = g = e and c = 1. Equilibrium E 5 The candidate Lyapunov function here is set as follows: Differentiating and using (10), we obtain Global stability now follows by setting f = g = e and c = 1. Equilibrium E 6 In this case the candidate is so that, using (12), We finally choose f = g = e and c = 1 to obtain negative definiteness. Equilibrium E 7 The candidate function is here The choice f = g = e and c = 1 renders A negative definite, and thus, global stability is achieved. Bifurcations Bifurcation E 0 − E 1 At E 0 , imposing the critical value r 0,1 = m, one eigenvalue vanishes and E 0 = E 1 . Now v = (1, 0, 0) T is the eigenvector of J (E 0 , r 0,1 ), A transcritical bifurcation then exists between E 0 and E 1 for r 0,1 = m in view of: Bifurcation E 0 − E 2 Here we impose the threshold value s 0,2 = n, to get a vanishing eigenvalue and E 0 = E 2 . The eigenvectors of J (E 0 , s 0,2 ) and its transpose are, respectively, v = (0, 1, 0) T and w = (0, 1, 0) T . A transcritical bifurcation occurs because of Bifurcation E 0 − E 3 In this case we set u 0,4 = p and find E 0 = E 4 . The required eigenvalues of J (E 0 , u = p) and [J (E 0 , u 0,4 )] T are v = (0, 0, 1) T and w = (0, 0, 1) T . The transcritical bifurcation arises in view of Bifurcation E 1 − E 3 To verify the transcritical bifurcation in this case to have E 1 = E 3 , we set The required eigenvectors are then v = (e(s − n), r − m, 0) T and w = (0, 1, 0) T , and the following conditions hold: The eigenvectors of the Jacobian and its transpose at this point with the parameter choice b = b 1,5 are, respectively, v = (e(u − p), 0, r − m) T and w = (0, 0, 1) T . The transcritical bifurcation is obtained in view of the conditions Bifurcation E 2 − E 3 Here to obtain E 2 = E 3 we choose a 2,3 = (m − r )(eS 2 ) −1 and find the eigenvectors v = (e(n − s), r − m, 0) T and w = (1, 0, 0) T , as well as the conditions F a (E 2 , a 2,3 ) = (0, 0, 0) T , w T F a (E 2 , a 2,3 ) = 0, w T (D F a v) = e 2 S 2 (n − s) = 0, This indicates once more that a transcritical bifurcation indeed arises between these points. Bifurcation E 2 − E 6 Considering u 2,6 = u 0,4 = p, we have E 2 = E 6 and the eigenvectors v = (0, 0, 1) T and w = (0, 0, 1) T , and the conditions F u (E 2 , u 2,6 ) = (0, 0, 0) T , w T F u (E 2 , u 2,6 ) = 0, w T (D F u v) = 1 = 0, w T (D 2 F(E 2 , u 2,6 )(v, v)) = −2c L L = 0, thus showing the occurrence of the transcritical bifurcation.