Keywords: Aggregation methods, Asymptotic series expansion, Canard, Geometrical Singular Perturbation theory, Predator-prey models Mathematics Subject Classification: 34E17 -37G15 -92D25 -92D40

Two predator-prey model formulations are studied: for the classical Rosenzweig-MacArthur (RM) model and the Mass Balance (MB) chemostat model. When the growth and loss rate of the predator is much smaller than that of the prey these models are slow-fast systems leading mathematically to singular perturbation problem.

In contradiction to the RM-model, the resource for the prey are modelled explicitly in the MB-model but this comes with additional parameters. These parameter values are chosen such that the two models become easy to compare. In both models a transcritical bifurcation, a threshold above which invasion of predator into prey-only system occurs, and the Hopf bifurcation where the interior equilibrium becomes unstable leading to a stable limit cycle. The fast-slow limit cycles are called relaxation oscillations which for increasing differences in time scales leads to the well known degenerated trajectories being concatenations of slow parts of the trajectory and fast parts of the trajectory. In the fast-slow version of the RM-model a canard explosion of the stable limit cycles occurs in the oscillatory region of the parameter space.

To our knowledge this type of dynamics has not been observed for the RM-model and not even for more complex ecosystem models. When a bifurcation parameter crosses the Hopf bifurcation point the amplitude of the emerging stable limit cycles increases. However, depending of the perturbation parameter the shape of this limit cycle changes abruptly from one consisting of two concatenated slow and fast episodes with small amplitude of the limit cycle, to a shape with large amplitude of which the shape is similar to the relaxation oscillation, the well known degenerated phase trajectories consisting of four episodes (concatenation of two slow and two fast).

1

The canard explosion point is accurately predicted by using an extended asymptotic expansion technique in the perturbation and bifurcation parameter simultaneously where the small amplitude stable limit cycles exist. The predicted dynamics of the MB-model is in a large part of the parameter space similar to that of the RM-model. However, the fast-slow version of MB-model does not predict a canard explosion phenomenon.

Modelling, singular perturbation and bifurcation analyses of bitrophic food chains B.W. Kooi * and J.C. Poggiale * * ,

Introduction

There is already a long tradition in modelling predator-prey systems under various environmental conditions. In some realistic cases there are differences in the order of magnitudes of the ingestion and growth rates of the two populations leading to so called fast-slow dynamical systems. Various mathematical analyses techniques have been used to analyse these models. Here we focus on the application of perturbation techniques when the processes of the two population occur at different time-scales.

The starting point is the Rosenzweig-MacArthur (RM) model [START_REF] Rosenzweig | Graphical representation and stability conditions of predator-prey interactions[END_REF] where the prey population grows logistically and the predator-prey interaction is described by the Holling type II functional response. This predator-prey system is mathematically described by two ordinary differential equations, (ode)s, for the prey and the predator population. This model does not obey mass conservation law, however. Therefore we study also a version where additionally the prey population consumes abiotic nutrients. The resulting model is called the Mass Balance (MB) model where mass conservation is obeyed. This type of model is common practice in modeling food chains in a chemostat reactor which have been studied extensively: we refer to [START_REF] Smith | The Theory of the Chemostat[END_REF][START_REF] Kooi | Consequences of population models on the dynamics of food chains[END_REF]. These systems consist minimally of a nutrient, a prey and a predator. Consequently, the description of the dynamics involves three odes.

The resulting three dimensional model is reduced to a two dimensional model for prey and predator populations when realistic assumptions are made.

For both types of models it is well known that the long-term dynamics is either a stable equilibrium (boundary: prey only, or interior: predator-prey) or a stable predator-prey limit cycle. For the RM-model, in [START_REF] Hsu | On global stability of a predator-prey system[END_REF] it was proved that the coexistence equilibrium is globally stable and in [START_REF] Cheng | Uniqueness of a limit cycle for predator-prey system[END_REF] that the periodic solution is unique thus a globally stable limit cycle and in [START_REF] Yu | Elements of Applied Bifurcation Theory[END_REF] a detailed bifurcation analysis is performed.

Nutrient enrichment leads to the risk of extinctions when fluctuating densities during the limit cycle reach low values, called the "paradox of enrichment" [START_REF] Rosenzweig | Paradox of enrichment: destabilization of exploitation ecosystems in ecological time[END_REF]. Mathematically the first transition occurs at a transcritical bifurcation and the second at a Hopf bifurcation. This holds for both types of model the RM and MB-model.

When the growth and loss rate of the predator is much smaller than that of the prey, the system becomes a so called slow-fast system. These rates for the predator are multiplied by a perturbation parameter. These systems have been studied intensively in the literature, the RM-model with fast-slow dynamics in [START_REF] Muratori | Low-and high-frequency oscillations in three-dimensional food chain systems[END_REF][START_REF] Rinaldi | Slow fast limit-cycles in predator prey models[END_REF][START_REF] Rinaldi | Destabilizing factors in slow-fast systems[END_REF][START_REF] Deng | Food chain chaos with canard explosion[END_REF][START_REF] Lobry | Effect of population size in a predator-prey model[END_REF]. Slow-fast systems can naturally be identified, are often encountered in ecology. This is almost a rule when behavioral and population dynamics are considered at the same time. But also interactions between two populations like plankton and fish, plants and insects, herbs and trees are classical examples of systems with two time scales (see [START_REF] Rinaldi | Destabilizing factors in slow-fast systems[END_REF]). Another example where there are differences in the order of magnitudes of the ingestion and growth rates is a food chain of sewage-bacterium-worms often found in wastewater treatment plants (see [START_REF] Ratsak | Modelling the growth of an oligocheate on activated sludge[END_REF]). The worm is the water nymph Nais elinguis, a oligochaete species. In a previous paper [START_REF] Kooi | Aggregation methods in food chains[END_REF] we toke advantage of these different time scales to apply aggregation methods in order to simplify the models for the dynamics of the wastewater treatment plant system.

Using a perturbation technique often the slow variables are frozen with the calculation of the equilibria of the fast system but better approximations are obtained by asymptotic expansions [START_REF] Hoppensteadt | Analysis and Simulation of Chaotic Systems[END_REF][START_REF] Hek | Geometric singular perturbation theory in biological practice[END_REF]. When the perturbation parameter becomes small the resulting limit cycle is called the relaxation oscillation, see [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF], where the periodic orbits are phase plane curves with both fast and slow parts of the trajectory.

When the time-scales differences are very large, that is with small positive perturbation parameter, the presence of the Hopf bifurcation indicates the possible occurrence of so called canards. In physics this has been studied extensively for the van der Pol equations [START_REF] Van Der Pol | On relaxation oscillations[END_REF][START_REF] Eckhaus | Relaxation oscillations including a standard chase on french ducks[END_REF][START_REF] Diener | The canard unchained or how fast/slow dynamical problems bifurcate[END_REF][START_REF] Canalis-Durand | Formal expansion of van der pol equation canard solutions are gevrey[END_REF][START_REF] Dumortier | Canard cycles and center manifolds[END_REF][START_REF] Krupa | Geometric analysis of the singularly perturbed fold[END_REF][START_REF] Krupa | Relaxation oscillation and canard explosion[END_REF][START_REF] Brons | An iterative method for the canard explosion in general planar systems[END_REF]. More recently, in [START_REF] Poggiale | Study of a virus-bacteria interaction model in a chemostat:application of geometrical singular perturbation theory[END_REF][START_REF] Hek | Geometric singular perturbation theory in biological practice[END_REF] geometric singular perturbation techniques [START_REF] Fenichel | Geometric singular perturbation theory[END_REF] have been studied for application in biological practice, including these slow-fast predator-prey systems. In this theory invariant manifolds play an important role in the study of structural stability of dynamical systems or, when a degeneracy occurs, in understanding the nature of bifurcations. A trajectory is an example of an invariant manifold. Using this approach we will find that similar to the van der Pol equations case also for the RM-model application of asymptotic expansions [START_REF] Canalis-Durand | Formal expansion of van der pol equation canard solutions are gevrey[END_REF][START_REF] Brons | An iterative method for the canard explosion in general planar systems[END_REF] for the critical manifolds leads to a divergent expansion. Nevertheless we will show that we are still able to compute for which bifurcation parameter values a canard explosion occurs.

Bifurcation analysis results show that in the MB-model the singularity in the limiting time-scale differences is completely different from that in the RM-model where the prey grows logistically. Only local bifurcations occur and therefore continuation of equilibrium and limit cycles gives a full picture of the long-term dynamics.

The RM bitrophic food chain model

A standard two-level food chain model from the field of theoretical biology is the scaled Rosenzweig-MacArthur system (1963). The RM-model reads in dimensionless form derived in the Appendix A:

dx 1 dt = x 1 (1 -x 1 - a 1 x 2 1 + b 1 x 1 ) , ( 1a
)
dx 2 dt = x 2 ( c 1 a 1 x 1 1 + b 1 x 1 -d 1 ) , (1b) with x i (t) ∈ R + , t ≥ 0, i = 1, 2
, respectively the size of the prey and predator population.

The first part of (1a) models that the prey population grows logistically in absence of the predator. The hyperbolic relationship a 1 x 1 /(1+b 1 x 1 ) is called the Holling type II functional response. This expression occurs in both equations and models the consumption of the prey population by the predator population. Parameter a 1 > 0 is the searching rate for the prey and b 1 > 0 the product of the handling time of the prey by the predator and the searching rate. The ratio of the first term of (1b) and the last term on the right-hand sides of (1a) is the efficiency c 1 > 0. The last term of (1b) is death rate d 1 > 0 of the predator population but it can also model other biological processes of which the loss rate is proportional to the population size similar to maintenance.

For a description of the state variables and the biological meaning of the parameters of the predator-prey model the reader is referred to Table 1.

Table 1: List of parameters and state variables and their reference values. We take the efficiency c 1 and death rate d 1 both proportional to the positive perturbation parameter ε. For numerical studies we take that parameter a 1 co-varies with b 1 via a 1 = 5/3 b 1 , hence handling time is 3/5 and the values for parameters c = d = 1 for the RM-model and D = e 1 = 1 for the MB-model. This is without loss of generality.

parameter ref.values Interpretation t [0, ∞) Fast time variable τ [0, ∞) Slow time variable x 0 [0, ∞) Nutrient density x 1 [0, ∞) Prey biomass density x 2 [0, ∞) Predator biomass density x r (0, ∞) Nutrient concentration in reservoir a 1 a 1 = 5/3 b 1 Searching rate b 1 3,4 or 8 Handling time × searching rate c 1 = εc c = 1 Conversion efficiency d 1 = εd d = 1 Death rate e 1 1 Conversion efficiency D 1 = εD D = 1 Dilution rate ε [0, 1] Perturbation parameter
The set of equations analysed extensively in the literature that form a model with slow-fast dynamics reads

dx 1 dt = f (x 1 , x 2 , ε) = x 1 (1 -x 1 - a 1 x 2 1 + b 1 x 1 ) , ( 2a
)
dx 2 dt = εg(x 1 , x 2 , ε) = εx 2 ( a 1 x 1 1 + b 1 x 1 -1) , (2b) with x i ∈ R, i = 1, 2.
We introduced c 1 = εc and d 1 = εd and toke for both parameters their reference value 1. The efficiency is again the ratio of the first term of (2b) and the last term on the right-hand sides of (2a) and is now equal to ε.

The functions f : R 3 → R and g : R 3 → R are of class smooth enough. The time-scale separation parameter ε is introduced in the model to implement trophic time diversification.

For ε 1 this is called a fast-slow system.

In the mathematical literature, factor ε is treated as a perturbation parameter, justified and described by the ratio between the linear death rate of the predator and the linear growth rate of the prey in [START_REF] Rinaldi | Slow fast limit-cycles in predator prey models[END_REF][START_REF] Hek | Geometric singular perturbation theory in biological practice[END_REF]. That is, only when the prey reproduce much faster than the predators and the predator is, in comparison, not so efficient, when the ratio ε becomes a small parameter. For, the efficiency c 1 = εc, where c is of order 1 (with a reference value 1) is proportional to the perturbation parameter ε.

The Holling type II hyperbolic relationship is derived by a time scale argument using a time budget modelling spend on searching for and handling of prey individuals by a predator individual. The ratio of these two terms is called the assimilation efficiency in ecology literature and the yield in the microbiology literature. When the units of both state variables equal then consequently ε < 1 means that there is a smaller than 100% biomass conversion, as is always the case in nature. We assume that the formed products during this conversion process, have no effects on other processes underlying model [START_REF] Brons | An iterative method for the canard explosion in general planar systems[END_REF]. In general, however, a very small conversion efficiency is not supported in the literature.

Observe that also the predator loss rate is multiplied by the same factor ε in order to facilitate coexistence. In other words when predators efficiency is low they also have to have a low loss rate in order to survive. Therefore the parameter ε affects two processes, mass conversion from prey biomass into predator biomass and the predator loss rate.

Existence and stability analysis of equilibria and limit cycles

In model ( 2) there are only three free parameters, namely a 1 , b 1 and ε which scales the efficiency and predator loss rate. The following stability analysis is classical for ε > 0 and therefore we recall some results regarding the dynamics of (2) and report some interesting results when ε → 0 and ε = 0

The one-parameter bifurcation diagram with varying b 1 where ε = 1 and a 1 co-varies with b 1 via a 1 = 5/3 b 1 , is shown in Fig. 1 for parameter b 1 as free parameter. The three relevant equilibria zero-E 0 , boundary-E 1 and interior equilibrium E 2 and the limit cycle L 2 are summarized in Table 2 and the bifurcation curves the transcritical bifurcation T C

and Hopf bifurcation H in Table 3. In [START_REF] Hsu | On global stability of a predator-prey system[END_REF] it was proved that the coexistence equilibrium E 2 is globally stable and in [START_REF] Cheng | Uniqueness of a limit cycle for predator-prey system[END_REF] that the periodic solution is unique thus a globally unique stable limit cycle L 2 . The numerical bifurcation results show that the limit cycle for parameter values above the supercritical Hopf bifurcation is stable and that the minimum values become very small for large b 1 . This phenomenon is related to the "paradox of enrichment" [START_REF] Rosenzweig | Paradox of enrichment: destabilization of exploitation ecosystems in ecological time[END_REF], because extinction due to stochastic fluctuations is likely.

For the stability analysis of equilibria we need the Jacobian matrix evaluated at point (x 1 , x 2 ) ( [START_REF] Kot | Elements of Mathematical Ecology[END_REF]), which reads

J = 1 -x 1 -a 1 x 2 1+b 1 x 1 + x 1 ( a 1 x 2 b 1 (1+b 1 x 1 ) 2 -1) -a 1 x 1 1+b 1 x 1 εx 2 ( a 1 1+b 1 x 1 -a 1 x 1 b 1 (1+b 1 x 1 ) 2 ) ε( a 1 x 1 1+b 1 x 1 -1) . ( 3
)
In E 0 the Jacobian matrix is

J 0 = 1 0 0 -ε . ( 4

Equilibria

System composition RM-model Bifs.

E 0 = (0, 0) Extinction E 1 = (1, 0) Prey-only E 2 = (x * 1 , x * 2 ) = ( 1 a 1 -b 1 , a 1 -b 1 -1 (a 1 -b 1 ) 2 ) Prey-predator MB-model E 0 = (0, 0) Extinction E 1 = (1, 0) Prey-only E 2 = (x * 1 , x * 2 ) = ( 1 a 1 -b 1 , a 1 -b 1 -1 ε(a 1 -b 1 ) 2 +a 1 -b 1 ) Prey-predator
a 1 (b 1 , ε) b 1 (a 1 = 5/3 b 1 ) Interpretation RM-model T C a 1 = b 1 + 1 1.5 E 1 → E 2 H a 1 = b 1 (b 1 +1) b 1 -1 4 E 2 → L 2 MB-model T C a 1 = b 1 + 1 1.5 E 1 → E 2 H a 1 = 2εb 2 1 +1+ √ 4b 2 1 ε(ε+1)+1 2ε(b 1 -1) b 1 = 4ε+ √ 16ε 2 +15ε (2ε) E 2 → L 2
The eigenvalues are the diagonal elements and the E 0 is always unstable. At the boundary equilibrium E 1 the Jacobian matrix reads

J 1 = -1 -a 1 -b 1 -1 1+b 1 0 ε( a 1 -b 1 -1 a 1
) .

(

The eigenvalues are again the diagonal elements and E 1 is stable when a 1b 1 -1 < 0 and unstable when a At the equilibrium E 2 the Jacobian matrix reduces to

J 2 = a 1 (b 1 -1)-b 1 (1+b 1 ) a 1 ((a 1 -b 1 ) -1 ε( a 1 -b 1 -1 a 1 ) 0 . ( 6
)
The complex pair of eigenvalues λ 1,2 of the Jacobian matrix evaluated at the equilibrium

E 2 read λ 1,2 (b 1 , ε) = μ(b 1 ) ± iω(b 1 , ε) . ( 7
)
We calculated the following real μ ∈ R and imaginary ω ∈ R parts as functions of parameter ε and b 1 where

a 1 = 5/3 b 1 Δ(b 1 , ε) = (Tr J 2 ) 2 -4 det J 2 = -144 + 72b 1 -9b 2 1 -60b 1 ε + 40εb 2 1 25b 2 1 , ( 8
)
μ(b 1 , ε) = 3(b 1 -4) 10b 1 ± 1 -sgn Δ(b 1 , ε) 2 |Δ(b 1 , ε)| , ( 9
) In Fig. 2 the curve where Δ(b 1 , ε) = 0 is plotted in the (b 1 , ε) parameter space. This curve separates regions where the equilibrium is a node Δ > 0 (two real eigenvalues) or a spiral or focus Δ < 0 (two conjugated complex eigenvalues). At Δ = 0 there are two equal real eigenvalues. From this figure we conclude that for ε = 1 and b 1 = 3 or b 1 = 8 the equilibria are foci while for ε = 0.01 the equilibria are nodes, stable and unstable, respectively.

ω(b 1 , ε) = ± 1 + sgn Δ(b 1 , ε) 2 |Δ(b 1 , ε)| . ( 10
)
The positions of the bifurcations T C and H are independent of the parameter ε. This is not true for the so called first Lyapunov coefficient 1 which determines whether the Hopf bifurcation is supercritical or subcritical, that is whether the emerging limit cycles are stable or unstable, respectively.

At the Hopf bifurcation with b 1 = 4 we have for ε > 0 a zero real part μ(4) = 0 and the positive imaginary part equal to the determinant of the Jacobian matrix J 2 : 3)] is satisfied. Using a similar procedure in Maple, [START_REF] Maple | Maple software. Maplesoft[END_REF] described in [START_REF] Yu | Elements of Applied Bifurcation Theory[END_REF], the first Lyapunov coefficient 1 is evaluated as

ω(4) = √ det J 2 = 1/2 √ ε,
1 = - 16 75 14ε + 1 ε 3/2 . ( 11
)
Since 1 is negative for ε > 0, the Hopf bifurcation is supercritical. Note that for lim ε →

0 we have 1 → -∞. That is, the first Lyapunov coefficient 1 becomes unbounded.
Furthermore, in the limiting case ε = 0 at the Hopf bifurcation where b 1 = 4 with μ(4) = 0

we get together with ω(4) = 1/2 √ ε = 0 that this point becomes a degenerated bifurcation where both eigenvalues are zero as in a Bogdanov-Takens bifurcation point. Exploration of these facts can be done with a blow-up technique [START_REF] Dumortier | Canard cycles and center manifolds[END_REF][START_REF] Dumortier | Geometric singular perturbation theory beyond normal hyperbolicity[END_REF] which is beyond the scope of this paper.

Phase-space analysis

In this section we discuss results in the phase-space of simulation in time. These results

where in all simulations the same initial conditions are used, are shown in Fig. 3 and Fig. 4:

• Left panels of Fig. 3: for b 1 = 3, in the region between the transcritical T C and Hopf H bifurcation in Fig. 1 where equilibrium E 2 is stable,

• Fig. 4: for b 1 = 4 at the Hopf bifurcation,

• Right panels of Fig. 3: for b 1 = 8 above the Hopf bifurcation in Fig. 1 where equilibrium E 2 is unstable.

In the top panels of Fig. 3 and Fig. 4 the f-nullcline where f (x 1 , x 2 , ε) = 0:

x 2 = (1 -x 1 )(1 + b 1 x 1 )
/a 1 and the g-nullcline where g(x 1 , x 2 , ε) = 0:

x 1 = 1/(a 1 -b 1 ), are
shown. Both are independent of ε. The graph of the f-nullcline is part of a parabola where it intersects the horizontal axis x 2 = 0 at x 1 = 1 and of the g-nullcline is just a vertical line through the equilibrium point E 2 .

Stable interior equilibrium

With b 1 = 3, see Fig. 3a with ε = 1 and Fig. 3c with ε = 0.01, the equilibrium E 2 is stable and there is convergence to the stable point.

For low ε = 0.01 value shown in Fig. 3c, initially after starting above the nullcline the solution goes rapidly to the vertical axis. There is almost no prey population and hence the predator population diminishes approximately exponential with rate d 1 = 1:

dx 2 dτ = -x 2 , ( 12
)
where differentiation is with respect to the slow time variable τ = εt. The solution crosses the critical point x 2 = 1/a 1 the intersection with the nullcline and leaves the vertical axis moving fast toward the stable part of the nullcline. Because this happens below the critical point (intersection of parabola and vertical axis) this phenomenon is called a delayed bifurcation which is explained in [START_REF] Rinaldi | Slow fast limit-cycles in predator prey models[END_REF]. Eventually the system converges to the stable equilibrium E 2 given in Table 2.

Hopf interior equilibrium

The equilibrium E 2 , see Table 2, at the Hopf bifurcation H coincides with the top of the f-nullcline parabola

(x 1 , x 2 ) T = b 1 -1 2b 1 , (b 1 + 1) 2 4a 1 b 1 . ( 13
)
At the Hopf point we have the equilibrium E 2

(x * 1 , x * 2 ) H = (x 1 , x 2 ) T = 3 8 , 15 64 . ( 14
)
The simulation results for system (2) are shown in Fig. 4. These results show that the solution finally converges very slowly to this, sometimes called a weakly attracting, equilibrium [START_REF] Guckenheimer | Normal Forms, Bifurcations and Finiteness Problems in Differential Equations[END_REF] at the top of the f-nullcline denoted by [START_REF] Ginoux | Flow curvature mathod applied to canard explosion[END_REF].

Unstable interior equilibrium

With b 1 = 8, see Fig. 3b,d equilibrium E 2 is unstable and there is no convergence to E 2 but to a stable limit cycle L 2 .

For positive ε = 0.01 in Fig. 3d initially the dynamics is similar to that for the stable case shown in Fig. 3c. The dynamics on the unique stable limit cycle L 2 consists of four

M 1 0 M 0 0 g(x 1 , x 2 , 0) f (x 1 , x 2 , 0) T • x 1 x 2 1 0.8 0.6 0.4 0.2 0 0.3 T C 0

Degenerate phase curves

With b 1 > 4 the limits of the limit cycles L 2 , the periodic relaxation oscillations, for lim ε → 0 are called degenerated phase curves shown in Fig. 5, see [START_REF] Guckenheimer | Normal Forms, Bifurcations and Finiteness Problems in Differential Equations[END_REF] It consists of four concatenated episodes (families of periodic sets). On two parts of the trajectory the dynamics is slow (single arrow): one along the vertical axis where the predator population decreases slowly and the other are along parts of the parabolic f-nullcline where the prey population decreases while the predator population predator population increases slowly.

On the other two parts of the trajectory the prey changes fast (double arrows) and the predator is constant: one from point T toward the vertical axis and one from this vertical axis to the point on the parabolic f-nullcline. This type of dynamics is discussed intensively in the ecological literature, we mention [START_REF] Rinaldi | Slow fast limit-cycles in predator prey models[END_REF][START_REF] Lobry | Effect of population size in a predator-prey model[END_REF] and references therein. These two episodes are connected by two fast episodes along the horizontal lines where the prey population size increases or decreases fast while the predator population is constant.

Starting slowly on the vertical axis above the f-nullcline, during the first episode, the solution orbit crosses the critical point T C and leaves this axis fast horizontally toward the parabola. Then it continues, during the second episode, slowly and during the third episode again along the parabola passing the top of the parabola point T given in [START_REF] Ginoux | Flow curvature mathod applied to canard explosion[END_REF].

Note that since b 1 > 4 the equilibrium E 2 indicated by the open circle in Fig. 5 on the f-nullcline, does not coincide with this non-hyperbolic point T . After crossing this point the solution orbit, during the forth episode, moves vastly horizontally toward the vertical axis again above the f-nullcline, and so on.

Note that the dynamics of the degenerate phase curves shown in Fig. 5 are indeed close to this for the RM-model in Fig. 3d where b 1 = 8 in both cases. However for the Hopf bifurcation b 1 = 4 case the resulting dynamics for small ε 1 does not look like that of Fig. 5 and more importantly also not for b 1 slightly above the Hopf bifurcation. Therefore we focus in the next sections on the dynamics for parameter values where the unstable equilibrium E 2 is just above the Hopf bifurcation value.

Singular perturbation problem

In this section the singular perturbation technique is used to analyse the singular perturbation problem where ε = 0 in the RM-model ( 2). We will also analyse the quasi-steady state solution or the relaxation oscillation in detail. The reader is referred to [START_REF] Hoppensteadt | Analysis and Simulation of Chaotic Systems[END_REF][START_REF] Jones | Geometric singular perturbation theory[END_REF][START_REF] Kevorkian | Multiple Scale and Singular Perturbation Methods[END_REF] for introductions into perturbation analysis.

Heuristic introduction

We start with a short overview of singular perturbation techniques. Singular perturbation theory deals with systems of the original form (2) where ε > 0. When ε 1 the system is a fast-slow system.

With ε = 0 we have the fast system also called called the layer system:

dx 1 dt = f (x 1 , x 2 (0), 0) = x 1 (1 -x 1 - a 1 x 2 (0) 1 + b 1 x 1 ) , ( 15a
)
dx 2 dt = 0 . ( 15b
)
The predator populations remains constant hence the trajectories are the horizontal lines in Fig. 5.

With a change of time-scale, where τ = εt, we call the resulting system the slow system:

ε dx 1 dτ = f (x 1 , x 2 , ε) , ( 16a
)
ε dx 2 dτ = εg(x 1 , x 2 , ε) . ( 16b
)
After substitution of ε = 0 we get

0 = f (x 1 , x 2 , 0) , (17a) dx 2 dτ = g(x 1 , x 2 , 0) = x 2 ( a 1 x 1 1 + b 1 x 1 -1) . (17b)
This differential algebraic equation (dae) is called the reduced system. The trajectory is either part of the vertical axis where x 2 = 0 or part on the parabola, the other f-nullcline of the original system (2) in Fig. 5. However, in the fast-slow context these are the sets of equilibria of the fast layer-system (15) and then they are called critical manifolds where

x 1 is acting as a parameter.

These heuristic results suggest the following approach for dealing with the two different time scales. The first step consists in setting ε = 0 which gives the set of fast equilibria of the fast system (15) yielding the algebraic equation (17a). This is the critical manifold which is the set of equilibria either at the vertical axis or the parabola where

f (x 1 , x 2 , 0) = 0.
From the results presented below in this section the parabolic critical manifold has two branches, one stable on the right-hand side (solid curve) and one unstable on the left-hand side (dashed curve) of point T in Fig. 5.

With good hypothesis (see below for the details), on the part of the parabola,

f (x 1 , x 2 , 0) = 0 is equivalent to x 1 = p(x 2 )
and we can substitute x 1 by p(x 2 ) in equation (17b). The result is the slow or reduced system

0 = f (p(x 2 ), x 2 , 0) , dx 2 dτ = g(p(x 2 ), x 2 ) . ( 18
)
The differentiation is with respect to τ and the smooth functions f and g are defined in [START_REF] Brons | An iterative method for the canard explosion in general planar systems[END_REF]. Then the slow system becomes

dx 2 dτ = x 2 (a 1 -b 1 )p(x 2 ) -1 b 1 p(x 2 ) + 1 , x 1 = p(x 2 ) = 1 2b 1 b 1 -1 ± (b 1 + 1) 2 -4a 1 b 1 x 2 . ( 19
)
An alternative method is to use instead of x 1 = p(x 2 ) of which the dynamics is described by the solution of ( 19), the inverse function x 2 = q(x 1 ) also derived from (17a). This gives the relationship

x 2 = q(x 1 ) = 1 a 1 (1 -x 1 )(1 + b 1 x 1 ) , ( 20
)
and using (17b) the differential equation

dx 2 dτ = g(x 1 , q(x 1 )) = dq dx 1 dx 1 dτ ,
and we get formally

dx 1 dτ = q(x 1 )(a 1 x 1 -(1 + b 1 x 1 )) (1 + b 1 x 1 ) dq/dx 1 , dq dx 1 = 1 a 1 b 1 (1 -2x 1 ) -1 . ( 21
)
Note that this expression is zero at point T given in [START_REF] Ginoux | Flow curvature mathod applied to canard explosion[END_REF]. Hence the denominator of model [START_REF] Jones | Geometric singular perturbation theory[END_REF] at that point T is zero. When the numerator is unequal zero, this means that the rate of change becomes unbounded at T which is a singular point of model [START_REF] Jones | Geometric singular perturbation theory[END_REF]. Only when in the special case b 1 is the Hopf bifurcation parameter value, point T is a limit point where the numerator is also zero. This is studied in the Appendix B.

Geometric singular perturbation techniques

We discuss the singular perturbation problem outlined in the previous section for the case where ε is not zero but small and positive: 0 < ε 1. Here we follow the geometric singular perturbation techniques.

Let us consider system (2) again. For ε = 0 the f-nullclines, the set

{(x 1 , x 2 )|f (x 1 , x 2 , 0) = 0, x 1 ≥ 0, x 2 ≥ 0} consist of two types of critical manifolds,see Fig. 5 M 0 0 = (x 1 , x 2 )|x 1 = 0, x 2 ≥ 0 (22a) M 0 1 = (x 1 , x 2 )|x 2 = 1 a 1 (1 -x 1 )(1 + b 1 x 1 ), x 1 , x 2 ≥ 0 . ( 22b
)
They form a set of equilibria of the fast system system (15). In the previous section we studied the dynamics for ε = 0 and now we will consider 0 < ε 1.

To that end, let us remind the statement of Fenichel's theorem. We consider differential system of the form:

dX dt = F (X, Y, ε) , ( 23a
) dY dt = εG(X, Y, ε) , ( 23b
)
dε dt = 0 , ( 23c
)
where F and G are sufficiently smooth. We assume that the set F (X, Y, 0) = 0 can locally be written as X = H(Y ), which defines a critical manifold. If, for all Y in a given compact set D, all the eigenvalues of DF DX (H(Y ), Y, 0) have a non-vanishing real part, then the critical manifold is said normally hyperbolic. In this case, there exists ε 0 and a map p defined on

D × [0, ε 0 [ such that: i) H(Y ) = p(Y, 0);
ii) the graph of p is invariant under the flow associated to the original differential system [START_REF] Kooi | Modelling the dynamics of traits involved in fighting-predators-prey system[END_REF];

iii) the graph of p is tangent to the central space associated to the lineralisation of the system at (H(Y ), Y, 0).

As a consequence, both critical manifolds M 0 0 and M 1 0 are normally hyperbolic and there exists ε 0 such that for 0 < ε < ε 0 , there are locally invariant manifolds M 0 ε and

M 1 ε except in the neighborhood of point T (x 1 , x 2 ) = (x 1 , x 2
) and in the neighborhood of the intersection between M 0 0 and M 1 0 on the vertical axis. Indeed, at those points, the derivative of the fast part vanishes, which contradicts the assumptions of the theorem statement.

Using its invariance, the perturbed manifold M 1 ε can be approximated by asymptotic expansions in ε. It can be described as a graph

{(x 1 , x 2 )|x 2 = q(x 1 , ε), x 1 ≥ 0, x 2 ≥ 0} . ( 24
)
This manifold is invariant when the following equality holds

dx 2 dt = dx 2 dx 1 dx 1 dt = dq(x 1 , ε) dx 1 dx 1 dt , (25)
yields with equation ( 23) and x 2 = q(x 1 , ε):

∂q(x 1 , ε) ∂x 1 dx 1 dt = εq(x 1 , ε) x 1 (a 1 -b 1 ) -1 1 + b 1 x 1 . ( 26
)
Then [START_REF] Kooi | Modelling the dynamics of traits involved in fighting-predators-prey system[END_REF] gives with x 2 = q(x 1 , ε) the invariance condition

∂q(x 1 , ε) ∂x 1 x 1 1 -x 1 - a 1 q(x 1 , ε) 1 + b 1 x 1 = εq(x 1 , ε) x 1 (a 1 -b 1 ) -1 1 + b 1 x 1 , (27)
or using 1

+ b 1 x 1 > 0 ∂q(x 1 , ε) ∂x 1 x 1 (1 -x 1 )(1 + b 1 x 1 ) -a 1 q(x 1 , ε) = εq(x 1 , ε) x 1 (a 1 -b 1 ) -1 . ( 28
)

Asymptotic expansion

The following asymptotic expansion in ε is introduced:

q(x 1 , ε) = q 0 (x 1 ) + εq 1 (x 1 ) + ε 2 q 2 (x 1 ) + . . . , (29)
hence

∂q ∂x 1 = dq 0 dx 1 + ε dq 1 dx 1 + ε 2 dq 2 dx 1 + . . . . ( 30
)
Substitution into [START_REF] Krupa | Geometric analysis of the singularly perturbed fold[END_REF] gives

( dq 0 dx 1 + ε dq 1 dx 1 + ε 2 dq 2 dx 1 + . . .)x 1 (1 -x 1 )(1 + b 1 x 1 ) -a 1 (q 0 (x 1 ) + εq 1 (x 1 ) + ε 2 q 2 (x 1 ) + . . .) = ε(q 0 (x 1 ) + εq 1 (x 1 ) + . . .) x 1 (a 1 -b 1 ) -1 . ( 31
)
Gathering orders of ε results for O(1) and assuming x 1 > 0 in:

q 0 (x 1 ) = (1 -x 1 )(1 + b 1 x 1 ) a 1 , dq 0 dx 1 = b 1 -1 -2x 1 b 1 a 1 . ( 32
)
At b 1 = 4 we have x 1 = (b 1 -1)/(2b 1 ) and hence dq 0 /dx 1 = 0.

For O(ε) and using an updated form of ( 28)

∂q(x 1 , ε) ∂x 1 x 1 a 1 (q 0 -q(x 1 , ε)) = εq(x 1 , ε) x 1 (a 1 -b 1 ) -1 , (33)
gives

q 1 (x 1 ) = q 0 (x 1 ) (x 1 (a 1 -b 1 ) -1) -a 1 x 1 dq 0 dx 1 . ( 34
)
At b 1 = 4 the numerator and denominator are both zero. we have x 1 = (b 1 -1)/(2b 1 ) and

dq 0 /dx 1 = 0 but also since it is a equilibrium x 1 (a 1 -b 1 ) -1 = 0.
For O(ε 2 ) in [START_REF] Marvá | Reduction of slow-fast periodic systems with applications to population dynamics models[END_REF] gives

q 2 (x 1 ) =q 1 (x 1 ) dq 1 dx 1 x 1 a 1 + x 1 (a 1 -b 1 ) -1 -a 1 x 1 dq 0 dx 1 . ( 35
)
At point T we have

x 1 = x 1 = (b 1 -1)/(2b 1 )
where dq 0 /dx 1 = 0. Consequently, at that point the denominator in the expression for q i (x 1 ), i > 0 is zero. Therefore the coefficients q i (x 1 ), i > 0 are unbounded when the numerator is not equal zero. Only when the parameter b 1 = 4 is at the Hopf bifurcation the numerator is zero and the coefficients q i (x 1 ) remain finite.

Asymptotic expansion in phase space

The expression for q 0 describes the critical manifold M 1 0 . This expression is the inverse (when it exists) of p(x 2 ) in [START_REF] Hsu | On global stability of a predator-prey system[END_REF]. The voluminous expressions for the higher order q i , i > 1 coefficients obtained by equating the O(ε i ) terms on the left-and right-hand side of the invariance condition [START_REF] Krupa | Geometric analysis of the singularly perturbed fold[END_REF], are not given here but are available using Maple, [START_REF] Maple | Maple software. Maplesoft[END_REF]. This yields the approximation of the perturbed slow manifold M 1 ε For ε = 0 the limit (16b) prescribes the singular slow flow on M 1 0 with x 2 = q 0 ( x 1 )

given by ( 32)

d x 1 dt = x 1 (1 -x 1 - a 1 q 0 ( x 1 ) 1 + b 1 x 1 ) . ( 36
)
For sufficiently small non-zero ε 1 the flow on the perturbed slow manifold M 1 ε can be approximated by inserting x 2 = q( x 1 , ε) with q( x 1 , ε) given by [START_REF] Krupa | Relaxation oscillation and canard explosion[END_REF]. In order to simulate the model we solve

d x 1 dt = x 1 1 -x 1 - a 1 q( x 1 , ε) 1 + b 1 x 1 , ( 37
)
with properly chosen the initial values.

The results are shown in Fig. 6 for b 1 = 3, where the equilibrium E 2 is stable. They show that in this case the solution of the original model on M 1 ε is already well approximated by the second order approximation when ε = 0.1.

At the Hopf bifurcation

In Fig. 7c the graph of the function q(x 1 , ε) is shown for b 1 = 4. While the solution converges to the neutral stable equilibrium E 2 where the rate of convergence becomes slower when approaching the non-hyperbolic point (see also Fig. 4b,c) the asymptotic expansion approximation explodes and there is discontinuity at

x 1 = x * 1 = x 1 : lim x 1 ↓x 1 q(x 1 , ε) = ∞ and lim x 1 ↑x 1 q(x 1 , ε) = -∞ . ( 38
)
In the next section we extend the asymptotic expansion in ε by varying parameter b 1 in addition to ε in order to repair this unwanted property.

Canard explosion

In this section we analyse the Canard dynamics that occurs for b 1 values just above the Hopf bifurcation at b 1 = 4 similar to the analysis performed in [START_REF] Brons | An iterative method for the canard explosion in general planar systems[END_REF]. Other papers on canards are [START_REF] Eckhaus | Relaxation oscillations including a standard chase on french ducks[END_REF][START_REF] Diener | The canard unchained or how fast/slow dynamical problems bifurcate[END_REF][START_REF] Dumortier | Canard cycles and center manifolds[END_REF][START_REF] Lobry | Effect of population size in a predator-prey model[END_REF][START_REF] Brons | An iterative method for the canard explosion in general planar systems[END_REF]. We will expand the asymptotic expansion discussed in the previous section where the equilibrium point is unstable and the system itself shows oscillatory behaviour, a stable limit cycle.

We start with an exploration of this oscillatory behaviour for b 1 values just above the Hopf bifurcation by simulation of the full model in time for small ε values. In order to study the dynamics for small ε in more detail we re-analyse the continuation as in Fig. 1 where ε = 1 in the region close to the Hopf bifurcation at b 1 = 4. The results are shown for ε = 0.01 in Fig. 8. The continuation of the curves calculated with AUTO [START_REF] Doedel | Auto 07p: Continuation and bifurcation software for ordinary differential equations[END_REF] failed for small ε. This is due to the part of the cycle close to the vertical axis where x 1 is small, see Fig. 3bd. In order to avoid this dynamics we study now an augmented system

dx 1 dt = δ + f (x 1 , x 2 , ε) = δ + x 1 (1 -x 1 - a 1 x 2 1 + b 1 x 1 ) , ( 39a
)
dx 2 dt = εg(x 1 , x 2 , ε) = εx 2 ( a 1 x 1 1 + b 1 x 1 -1) , ( 39b
)
where δ is a small allochthonous input rate of the prey population. Addition of this extra term removes the transcritical bifurcation at x 2 = 1/a 1 because it is structurally unstable with respect to such a perturbation. show that the canard is a robust and smooth phenomenon that occurs for stable limit cycles. When ε ↓ 0 the canard explosion point converges to the point T .

In Fig. 11 the shape of the limit cycles is show for two values of b 1 just below and above the sudden changes: b 1 = 4.04019 and b 1 = 4.04061 where ε = 0.01. These results show how the unique stable limit cycle changes shape very abruptly at b 1 ≈ 4.0403 where b 1 is varied keeping ε fixed: the canard explosion. In Fig. 12 the shape of the limit cycles is show for four values: ε = 0.008, 0.009, big stable limit cycles and ε = 0.01, 0.011, small stable limit cycles with fixed b 1 = 4.0403, the value calculated with the extended asymptotic expansion technique. These results show the dynamics with ε values above the canard explosion at approximately ε = 0.01, are limit cycles with small amplitudes consisting of two concatenated slow (close to the critical manifold and just above the parabola) and fast (almost horizontal part of the trajectory below the parabola). Below the canard explosion the limit cycles with large amplitudes consist of four episodes, two slow (close to the critical manifold, one just above the parabola and one close to the vertical axis) and two fast (almost horizontal, one leaving close from point T and one leaving the vertical axis and landing close to the parabola).

Carefully examination indicates that there is one point where the trajectories for the four ε values intersect. Starting at that point and changing ε gives the unique stable limit 39) with prey input rate δ = 0.0001. The parameter continuation is here successful to very small values of ε.

cycles calculated and this shows that these curves change smoothly, only the sensitivity of the shape of the cycles is very large at the explosion point.

Asymptotic expansion and canard explosion

Following [START_REF] Ginoux | Flow curvature mathod applied to canard explosion[END_REF] we repeat the extended asymptotic expansion in ε technique introduced in (29) now near the Hopf bifurcation point. To that end the expansion is not only taken for the function x 2 = q(x 1 , ε) evaluated at q 0 (x 1 ) (see [START_REF] Maple | Maple software. Maplesoft[END_REF]) but also in the bifurcation parameter b 1 = b 1 (ε) and consequently a 1 (ε) = 5/3 b 1 (ε) evaluated at the Hopf bifurcation point b 10 = 4 with a 01 = 20/3. This comes with more freedom which leads to an extra criterion: the singularity of the approximation at point T is removed. This function is now denoted as r(x 1 , ε) and the expansion is evaluated at r(

x 1 = x 1 , ε = 0). When b 1 = b 10 we have r(x 1 , ε) = q(x 1 , ε).
Similar to the invariance condition of the perturbed slow manifold M 1 ε we derive now the adapted equation that replaces [START_REF] Krupa | Geometric analysis of the singularly perturbed fold[END_REF] where we substitute a 1 (ε) = 5/3 b 1 (ε):

∂r ∂x 1 x 1 (1 -x 1 )(1 + b 1 (ε)x 1 ) -5/3 b 1 (ε)r(x 1 , ε) = εr(x 1 , ε) 2/3 x 1 b 1 (ε) -1 . ( 40
)
The following extended asymptotic expansion in ε for r(x 1 , ε) is introduced as follows: and

x 2 = r(x 1 , ε) = r 0 (x 1 ) + εr 1 (x 1 ) + ε 2 r 2 (x 1 ) + . . . , (41)
∂r ∂x 1 = dr 0 dx 1 + ε dr 1 dx 1 + ε 2 dr 2 dx 1 + . . . , ( 42
)
whereby b 1 (ε) is given by b 1 (ε) = b 10 + εb 11 + ε 2 b 12 + ε 3 b 13 + ε 4 b 14 . . . , ( 43
)
where r j and b 1j , j = 1 • • • are independent of ε and are described by the invariance condition ( 40) by equality order by order of powers of ε of this condition. Equating O [START_REF] Auger | Aggregation methods in dynamical systems and applications in population and community dynamics[END_REF] terms yields:

r 0 = (1 -x 1 )(1 + b 10 x 1 ) 5/3 b 10 . ( 44
)
With b 10 = 4 this gives the zero order approximation of M 1 ε equal to M 1 0 . Using [START_REF] Rosenzweig | Graphical representation and stability conditions of predator-prey interactions[END_REF] and substitution of ( 41), ( 42) and ( 43) in [START_REF] Ratsak | Modelling the growth of an oligocheate on activated sludge[END_REF] and equating the resulting first order O(ε) terms, yields: From this point a recursive procedure can be followed. One by one, r i , i ≥ 2 is determined by taking the i th order term in the invariance condition [START_REF] Ratsak | Modelling the growth of an oligocheate on activated sludge[END_REF] equal to zero and thereafter the term b th 1(i-1) by the condition that this term is continuous at x 1 = x 1 = (b 10 -1)/(2b 10 ). This means that the free parameter is chosen such that the singularity is removed. The requirement that the sum of second order terms is zero gives

r 1 (x 1 ) = 3(1 -x 1 )(1 -b 10 + 2x 1 b 10 )b 11 x 1 b 10 -x 1 b 2 10 -3b 10 + 2x 2 1 b 3 10 ) 5b 2 10 (1 + 2x 1 b 10 -b 10 )x 1 . ( 45
r 2 = (1 -x 1 ) (-288x 4 1 + 108x 3 1 )b 12 + (72x 4 1 -27x 3 1 )b 2 11 + (96x 3 1 + 84x 2 1 )b 11 + 256x 3 1 -128x 2 1 -112x 1 -16 960x 3 1 (8x 1 -3) . ( 46
)
At the support point T being the Hopf bifurcation point with b 10 = 4 and x 1 = x 1 given by ( 14) we have (-288x 4 1 + 108x 3 1 )b 12 = 0. Furthermore the denominator is zero b 10 = 4. The expression for b 11 is then obtained by setting also the numerator equal to zero in order to remove the singularity at that support point. This is related to the case analysed in the previous section for the asymptotic expansion

x 2 = q(x 1 , ε) at q(x 1 , ε) for coefficient q 1 given in [START_REF] Monod | Recherches sur la croissance des cultures bactériennes[END_REF] where this expression was also unbounded. Now, because we have more freedom we can take b 11 so that the expression stays bounded at this point. The expression b 11 = 100/27 is obtained by substitution of the equilibrium value for x 1 = x 1 = 3/8 given by ( 14) into [START_REF] Smith | The Theory of the Chemostat[END_REF].

Further recursion gives higher order approximations, again first r n with b n-1 then r n+1

with b n and so on. We calculated the following fourth order approximation

b 1 (ε) = b 10 + εb 11 + ε 2 b 12 + . . . , b 1 (ε) = 4 + ε 100 27 + ε 2 58700 2187 + ε 3 80536900 177147 + ε 4 171270040300 14348907 . ( 47
)
In Fig. 11 besides the shape of the limit cycles for two values b 1 = 4.04019 (small limit cycle) and b 1 = 4.04061 (big limit cycle) also the result of the extended asymptotic expansion r(x 1 , ε), where ε = 0.01, (41) with b 1 = 4.0403 is shown. The approximation for M 1 ε follows the limit cycle closely up-to a separation point where the small limit cycle bends toward the nullcline which is crossed where the rate of x 1 changes sign and the rate becomes fast. This occurs for the smallest b 1 = 4.04019-value. For the big limit cycle at b 1 = 4.04061, from the separation point the trajectory along the cycle continuous to move toward the vertical axis. The asymptotic expansion for M 1 ε with b 1 ≈ 4.0403 continuous after the separation point between the two limit cycles before the approximation becomes unbounded.

These results show how the unique stable limit cycle changes its shape very abruptly when b 1 ≈ 4.0403 is varied keeping ε fixed, here in our example 0.01: the canard explosion. That the iteration process converges to this bifurcation parameter value where the approximation of the asymptotic expansion r(x 1 , ε) works for the limit cycle with the small amplitude that intersect with the parabolic f-nullcline vertically at the minimum predator size during the limit cycle. This makes in plausible that the iteration procedure yields indeed the canard point.

However, the extended asymptotic expansion r(x 1 , ε) (41) is divergent as shown in Fig. 13 where the coefficients b 1i as function of i is depicted.

It was shown in [START_REF] Ramis | Poincaré et les développements asymptotiques[END_REF][START_REF] Ramis | Les développements asymptotiques après poincaré : continuit et[END_REF] and reference therein, that the summation up to the smallest term gives an optimal (and very accurate) approximation in the case of the van der Pol system. In Fig. 14, again with allochthonous prey input where δ = 0.0001, the parameter value where the explosion occurs is plotted. This b 1 parameter value is taken from Fig. 10.

Fortunately, the result presented are in agreement.

The MB bitrophic food chain model

This section presents the Monod chemostat model [START_REF] Monod | Recherches sur la croissance des cultures bactériennes[END_REF][START_REF] Smith | The Theory of the Chemostat[END_REF][START_REF] Kooi | Consequences of population models on the dynamics of food chains[END_REF][START_REF] Kooi | Aggregation methods in food chains with nutrient recycling[END_REF]. In this model the nutrients consumed by the prey are modelled explicitly instead of using a logistic growth model for the growth of the prey in absence of the predator. Let x 0 (t) denote the density of the nutrient, and x i (t) ∈ R + , t ≥ 0, i = 1, 2 the biomass densities of prey and predator, respectively. The scaled version of the Monod model reads

dx 0 dt = (x r -x 0 )D 1 -a 0 x 0 x 1 , ( 48a
)
dx 1 dt = a 0 x 0 x 1 -D 1 x 1 - a 1 x 1 x 2 1 + b 1 x 1 , ( 48b
)
dx 2 dt = e 1 a 1 x 1 x 2 1 + b 1 x 1 -D 1 x 2 . ( 48c
)
where x r is the concentration of nutrient in the reservoir and D 1 the dilution rate, the rate at which the nutrient enters and all trophic levels are exported from the chemostat reactor. The second term of (48a) and the first term of (48b) model the Lotka-Volterra functional response at the prey-nutrient level. The first term of (48c) and the last term on the right-hand sides of (48b) represent the Holling type II functional response. The efficiency e 1 is the ratio of these two terms. For simplicity we assume e 1 = 1 and a 0 = 1.

It can be shown that all solutions system (48) starting in the non-negative cone eventually lie in the set

Ω = x = (x 0 , x 1 , x 2 ) ∈ R 3 + : x 0 + x 1 + x 2 ≤ x r . ( 49
)
So, the asymptotic behavior of system (48) is bounded by x r . It is possible to decouple this system by the introduction of the function

H(t) = x 0 (t) + x 1 (t) + x 2 (t) -x r , t ≥ 0 . ( 50
)
This gives

dH dt = -D 1 H , ( 51a
)
dx 1 dt = H + x r -x 1 -x 2 x 1 -D 1 x 1 + D 1 a 1 x 1 x 2 1 + b 1 x 1 , ( 51b
)
dx 2 dt = D 1 x 2 a 1 x 1 1 + b 1 x 1 -1 . ( 51c
)
When furthermore H = 0 for the asymptotic dynamics we can study the two dimensional predator-prey system

dx 1 dt = x 1 x r -x 1 -x 2 -D 1 -D 1 a 1 x 2 1 + b 1 x 1 , ( 52
)
dx 2 dt = D 1 x 2 a 1 x 1 1 + b 1 x 1 -1 . ( 53
)
In order to be able to make a clear comparison with the RM-model formulations possible we use x r = 1 + D 1 . This gives with

D 1 = ε dx 1 dt = x 1 1 -x 1 -x 2 -ε a 1 x 2 1 + b 1 x 1 , ( 54a
)
dx 2 dt = εx 2 a 1 x 1 1 + b 1 x 1 -1 . ( 54b
)
We call this the MB-model. One main difference with the RM-model is that the last term of (54a) is proportional to ε and consequently the efficiency, is constant. Another difference with the RM-model is that the logistic prey growth equation x 1 (1x 1 ) is replaced by the

expression x 1 (1 -x 1 -x 2 )
with one extra term namely that of the predator biomass, x 2 .

In absence of the predator, x 2 = 0, both expressions for the prey growth are the same.

For the three trophic system including the nutrients, which is used with the derivation of the MB-model, the biomass allocated in the predator gives a feed-back mechanism so that there is less nutrient available for the prey. In the food chain the predator has two adverse effects on the growth of prey population. Firstly the prey is consumed by the predator and they consume building-block material not only for themselves but also for their predators population that can only exists when the prey exists in the absence of inter-guild predation.

The biological interpretation of the -D 1 H = -εH term in (51a) is the difference between the influx rate and the out-flux rate of the total biomass expressed in the biomass of the predator. In [START_REF] Smith | The Theory of the Chemostat[END_REF] (50) is used to show that Monod's model is dissipative and that the system converges asymptotically to the manifold H = 0 where the influx rate and out-flux rate of the total biomass are the same. Observe that where H = 0 there is with respect to the RM-model a new x 2 term in the expression for the prey growth rate.

Existence and stability analysis of equilibria and limit cycles

For the analysis of model is the chemostat environment we refer to [START_REF] Smith | The Theory of the Chemostat[END_REF]. The three relevant equilibria are summarized in Table 2. In Table 3 the bifurcation analysis results are given.

Important difference of these results with those for the RM-model is that while the expressions for the T C bifurcation are the same those for the Hopf H bifurcation still depends on ε. There is convergence to a stable equilibrium, similar as we found for the RMmodel in Fig. 3a,c.

Phase-space analysis

For b 1 = 8 the results are shown in Fig. 17. These results differ much from the RMmodel in Fig. 3b,d case where the equilibrium E 2 was unstable. Here this holds true for ε = 1 (Fig. 17a) but for smaller values the equilibrium becomes stable (Fig. 17b,c).

The degenerate phase point

With ε = 0 substituted in the MB-model (54) there is no input of nutrients (54a) and also no export of the abiotic and biotic elements from the reactor environment. The prey grows logistically to the equilibrium x 0 (0) + x 1 (0) and the predator population remains constant x 2 (0). Hence, the equilibrium E 2 is neutral stable. The degenerate phase curve is just this point which is an equilibrium point x 1 together with the initial predator size x 2 (0). This degenerate phase curve differs completely from that of the RM-model. This is a consequence of the fact that in the MB-model the second term of (54a) is proportional to ε and therefore the ratio of this and the first term of (54b), the efficiency, is constant.

Bifurcation analysis of MB-model

In order to find-out why this happens we calculated a two-parameter bifurcation diagram shown in Fig. 18 where besides b 1 (whereby a 1 = 5/3 b 1 ), parameter ε is the second variable. In this diagram the transcritical bifurcation T C, and the Hopf bifurcation H, curves are drawn (see Table 3 for the expressions that describe the curves).

The transcritical bifurcation T C, is the same for all models. This is obviously due to the fact that the model for the dynamics of the predator is the same for all models and furthermore that for the prey-only (x 2 = 0) equilibrium E 1 is also the same. For the MBmodel the Hopf bifurcation H terminates at the origin where b 1 → ∞ and ε → 0. There is a stable equilibrium E 2 in almost the whole b 1 > 4 range up to lim b 1 ↑ ∞ while in the RM-model there is a stable limit cycle L 2 .

From this we conclude that in the case of the MB-model the parameter ε can not to be used as a single perturbation parameter. We conclude that the complete model has to be analysed using a straight-forward phase-space and bifurcation analysis of the local bifurcations H and T C.

Discussion and conclusions

The use of time-scale separation technique has a long tradition in ecology and biochemistry, starting with the quasi-steady-state approximation (QSSA) used to derive the Holling types functional response [START_REF] Holling | Some characteristics of simple types of predation and parasitism[END_REF] and Michaelis-Menten kinetics.

In this paper we compared two fast-slow versions of predator-prey models: the RM- model and MB-model. In the classical RM-model the small perturbation parameter is proportional to the efficiency of the predator-prey trophic interaction and the predator death rate. In a series of papers [START_REF] Muratori | Low-and high-frequency oscillations in three-dimensional food chain systems[END_REF][START_REF] Rinaldi | Slow fast limit-cycles in predator prey models[END_REF][START_REF] Rinaldi | Destabilizing factors in slow-fast systems[END_REF] the dynamics of this bitrophic systems has been studied intensively.

In [START_REF] Lobry | Effect of population size in a predator-prey model[END_REF] this subject was also studied focusing on the dynamics close to the critical manifold M 0 0 (Fig. 5) (part of the vertical axis where the prey is absent) and related to a delayed (transcritical) bifurcation phenomenon. Recently in [START_REF] Hek | Geometric singular perturbation theory in biological practice[END_REF] the analysis technique based on the geometric singular perturbation theory was applied. This theory can also be used where the interior equilibrium is unstable and relaxation dynamics occurs.

We used the invariant manifold criterion together with an extended asymptotic expansion with respect to both, the perturbation parameter and the free bifurcation parameter. This gives an iteration process that approximates the M 1 ε slow manifold close to the parabolic critical manifold M 1 0 in the neighbourhood of its top. In the parameter range just above the Hopf bifurcation where stable limit cycles exist, this process converges to the point where the canard explosion occurs. In [START_REF] Brons | An iterative method for the canard explosion in general planar systems[END_REF] a similar method has been used: here the advantage is that all calculations are done in Maple, [START_REF] Maple | Maple software. Maplesoft[END_REF] with rational numbers.

Direct application, however, showed that the delayed transcritical bifurcation dynamics close to the critical manifold M 0 0 leads to an disturbing effect because the prey population 3. The transcritical bifurcation curve T C is for both models RM-model and MB-model the same. For ε = 1 the point of the Hopf bifurcations differ slightly but for lim ε ↓ 1 they differ essentially. In the MB-model there is a stable equilibrium E 2 in almost the whole b 1 > 4 range up to lim b 1 ↑ ∞ while there is a stable limit cycle L 2 in the RM-model. becomes very low. Therefore we introduced a small allochthonous input of prey in which case the non-generic transcritical bifurcation disappears. In [START_REF] Sugie | Uniqueness of limit cycles in a Rosenzweig-MacArthur model with prey immigration[END_REF] the uniqueness of the limit cycles in the RM-model with prey immigration is shown that supports the applicability of this additionally introduced mechanism.

The canard phenomenon found can be described as follows. Also for small perturbation parameter values the trajectories follow the stable limit cycles like in the case where no time scale differences occur ε = 1. However, just above the Hopf bifurcation point depending also on the perturbation parameter the stable limit cycle with small amplitude consists of two concatenated slow and fast episodes. In the phase space, the top part is just above the parabolic stable and unstable critical manifold and the bottom just below the horizontal line connecting the two points where the limit cycle intersects with the fnullcline vertically. At the canard point of the bifurcation parameter, the dynamics tends abruptly toward the relaxation dynamics. This transition point resembles what happens due to the delayed bifurcation effect where the trajectory leaves the vertical axis below the transcritical bifurcation T C in Fig. 5. For higher bifurcation parameters values the limit cycle with large amplitude changes smoothly and approaches the degenerated phase curves consisting of a concatenation of two slow and two fast episodes (Fig. 5).

Despite the approximate series expansion diverges, we found accurate approximations for small ε of the part the limit cycles originating from the Hopf bifurcation point. Furthermore, the numerical approximate asymptotic iterative scheme, gives very good approximations of the bifurcation parameter b 1 and perturbation parameter ε values where a canard explosion occurs (see Fig. 11). The non-dimensional model then reads

dx 1 dt = x 1 (1 -x 1 ) - a 1 x 1 1 + b 1 x 1 x 2 , ( 56a
)
dx 2 dt = x 2 ( c 1 a 1 x 1 1 + b 1 x 1 -d 1 ) , (56b)

B The slow dynamics on the parabola part of the fnullcline

We give the results for the slow dynamics in the degenerated case ε = 0 where the dynamics is on a part of the f-nullclines critical manifolds.

The slow dynamics on the vertical axis is described again by the system [START_REF] Fenichel | Geometric singular perturbation theory[END_REF]. We will now focus on the slow dynamics on the parabola shown in Fig. 19 where the dynamics of the reduced system ( 21) is solved numerically, yielding the prey population size x 1 (τ ) while ( 20) is used to get the associated predator population size x 2 (τ ). Figure 19a gives the results for the stable equilibrium E 2 , b 1 = 3 where a 1 = 5/3 b 1 = 5 case with initial condition x 2 (0) = 0.225. This value is below point T , where x 2 = 0.25. In this case there are two valid initial points given by [START_REF] Hsu | On global stability of a predator-prey system[END_REF]. Starting at the largest x 1 prey value on the right-hand parabolic branch there is convergence toward the equilibrium E 2 . But at the lowest prey value on the left-hand branch there is convergence to the trivial zero-solution where x 1 = 0 and x 2 = 0. Hence, the right-hand branch of the parabola is in the basin of attraction of the stable equilibrium of the reduced system equal to that of the original system E 2 . On the other hand for the left-hand branch it is in the basin of attraction of the zero solution. the equilibrium E 2 of the original system is unstable. In Fig. 19c starting at the lowest prey value and x 2 (0) = 0.15, the results are similar to that in the above discussed case. However, starting with the largest prey value on the critical manifold there is convergence to the limit point T for the reduced system, and not to the unstable equilibrium E 2 of the full model. Note that the vector field is not defined at the top T . In Fig. 19d the initial value is x 2 (0) = 0.16 where both simulations terminate at a limit point T of the reduced system. Hence, the unstable equilibrium E 2 of the full system is a separatrix between the two equilibrium points of the reduced system, being limit point T in [START_REF] Ginoux | Flow curvature mathod applied to canard explosion[END_REF] and the zero point E 0 (x 1 , x 2 ) = (0, 0).

Figure 1 :

 1 Figure 1: One-parameter bifurcation diagram for b 1 and a 1 = 5/3 b 1 of the RM-system (1) with ε = 1. The equilibria E 1 (below T C) and E 2 (between T C and H) as well as the maximum and minimum values for the limit cycle L 2 (above H) are shown. Note that the critical values where the bifurcations T C and H occur in this diagram are independent of ε since the expression for the real parts of the eigenvalues of the Jacobian matrix do not depend on ε.

1 -b 1 - 1 <

 111 0. The equality gives the transcritical bifurcation T C where b 1 = a 1 -1 or for the reference value a 1 = 5/3 b 1 we have b 1 = 3/2.

Figure 2 :

 2 Figure 2: Focus bifurcation where Δ = 0 curve for ε vs b 1 .

  for the two conjugated complex eigenvalues. The square root of the ratio μ(b 1 )/ω(b 1 ) measures the amplitude of the limit cycle that emerge for the Hopf bifurcation for b 1 > 4. Evaluation of this ratio just above b 1 = 4 gives that the oscillatory dynamics for values is as with the Hopf bifurcation till the canard explosion occurs. The derivative dμ/db 1 (b 1 )| b 1 =4 = 3/40 > 0 at the Hopf bifurcation equilibrium is positive and therefore the transversility condition for applying the normal form theorem [31, Theorem (3.

Figure 3 :Figure 4 :

 34 Figure 3: Phase-space analysis for system (2) describing the RM-model with b 1 = 3 (left panels) and b 1 = 8 (right panels) while a 1 = 5/3 b 1 . Top panels (a,b): ε = 1 and bottom panels (c,d): ε = 0.01. (a): Spiral stable equilibrium E 2 , (b): stable limit cycle L 2 , (c): node stable equilibrium E 2 and (d): stable limit cycle L 2 .

Figure 5 :

 5 Figure 5: Fast (double arrow) and slow (single arrow arrow) dynamics for system (2) describing the RM-model with b 1 = 8 where for lim ε → 0 the concatenated trajectories are the degenerated phase curves. The f-nullcline (parabola) and g-nullcline (vertical line through equilibrium) are shown. Point T is top of parabola given in (13) and T C is intersection point of f-nullcline and vertical axis.

Figure 6 :Figure 7 : 1 ε

 671 Figure 6: Results for the original system (2) (solid) describing the RM-model (29) with b 1 = 3 where a 1 = 5/3 b 1 and initial conditions x 1 = 0.86038 and x 2 (0) = 0.102646.

Figure 8 :

 8 Figure 8: One-parameter bifurcation diagram for b 1 and a 1 = 5/3 b 1 of the RM-system (1) with ε = 0.01.

Fig. 10 Figure 9 :

 109 Figure 9: One-parameter bifurcation diagram for ε with various b 1 = 4.01, 4.02, 4.03, 4.04, 4.0401, 4.0407, 4.05 values where a 1 = 5/3 b 1 of the RM-system (1). For b 1 > 4 the minimum and minimum populations values during the stable limit cycle ate show. Also the two dashed curves are shown for b 1 = 4.0409 and 4.04061. Between these two values the explosion occurs at ε = 0.01 see also Fig. 8. The curves terminates for low ε values. This is due to numerical problems of the parameter continuation for these low values. Theory predicts that these curves of maximums and minimums are continuing almost horizontally.

Figure 10 :

 10 Figure 10: One-parameter bifurcation diagram for ε with various b 1 = 4.00125, 4.0025, 4.005, 4.01, 4.02, 4.03, 4.04 values where a 1 = 5/3 b 1 of system (39) with prey input rate δ = 0.0001. The parameter continuation is here successful to very small values of ε.

Figure 11 :Figure 12 :

 1112 Figure 11: Phase-space diagram with ε = 0.01. Solid lines: stable limit cycle for two values of b 1 = 4.04019 (small cycle) and b 1 = 4.04061 (big cycle) while a 1 = 5/3 b 1 of the RM-system (1). Dashed line: extended asymptotic expansion r(x 1 , ε) of M 1 ε (41) where b 1 = 4.0403.

)

  It appears that the terms with b 11 , namely 3x 1 b 11 (b 10 -1 -2x 1 b 10 ) = 0, are already eliminated when they are evaluated at the Hopf bifurcation support point with b 10 = 4 and a 10 = 20/3 at the equilibrium E 2 with x 1 = x 1 = x * 1 = 1/(a 10b 10 ). At that point both the numerator and the denominator are zero but lim x 1 ↓x 1 r 1 (x 1 ) = lim x 1 ↑x 11 r 1 (x 1 and this shows that function r 1 (x 1 ) is continuous for all b 11 .

Figure 13 :Figure 14 :

 1314 Figure 13: Coefficients b 1i as function of i given by (47) with ε = 0.01.

Firstly we calculate by continuation of the parameter b 1

 1 the bifurcation diagram shown in Fig.15. This diagram looks very much the same as Fig.1for the RM model. The main difference is that the Hopf bifurcation occurs at a somewhat higher b 1 value.

Figure 15 :

 15 Figure 15: One-parameter bifurcation diagram for b 1 and a 1 = 5/3 b 1 of the MB-system (54) with ε = 1.

Figure 16

 16 Figure 16 displays the simulation results for the MB-model for three values of ε where b 1 = 3.There is convergence to a stable equilibrium, similar as we found for the RM-

Figure 16 :

 16 Figure 16: Phase-space analysis for system (54) describing the MB-model with a 1 = 5/3 b 1 , b 1 = 3, for three different values of ε: ε = 1, ε = 0.1 and ε = 0.01.

Figure 17 :

 17 Figure 17: Phase-space analysis for system (54) describing the MB-model with a 1 = 5/3 b 1 , and b 1 = 8, that is at the Hopf bifurcation point for three different values of ε: ε = 1, ε = 0.1 and ε = 0.01.

Figure 18 :

 18 Figure 18: Two-parameter bifurcation diagram with ε and 1/b 1 as free parameters. The expressions are given in Table3. The transcritical bifurcation curve T C is for both models RM-model and MB-model the same. For ε = 1 the point of the Hopf bifurcations differ slightly but for lim ε ↓ 1 they differ essentially. In the MB-model there is a stable equilibrium E 2 in almost the whole b 1 > 4 range up to lim b 1 ↑ ∞ while there is a stable limit cycle L 2 in the RM-model.

where a 1 =

 1 AK/R, b 1 = kAK, c 1 = C and d 1 = M/R.

Figure 19 :

 19 Figure 19: Phase-space analysis for slow system (19) of the RM-model. (a): At the top-panel left b 1 = 3 and a 1 = 5/3 b 1 , with initial values x 2 = 0.225 and x 2 = 0.256. The equilibrium E 2 is stable and the right-branch of the parabola is the basin of attraction. On the other hand the left-branch of the parabola is in the attraction basin of the zero solution. (b): At the top-panel right b 1 = 4 and a 1 = 5/3 b 1 , with initial values x 2 = 0.225. The equilibrium E 2 is unstable. (c) and (d): The lower panel with b 1 = 8 and a 1 = 5/3 b 1 the two initial values x 2 = 0.15 and x 2 = 0.16. The associated initial value of the x 1 (0) is for x 2 = 0.15 below the unstable equilibrium E 2 value. For x 2 = 0.16 the associated initial value of the x 1 (0) is above the unstable equilibrium E 2 value.

Figure

  Figure 19c,d were calculated with parameter values b 1 = 8 and a 1 = 5/3 b 1 = 40/3 where

Table 2 :

 2 Equilibria of RM-model (2) and MB-model (54).

Table 3 :

 3 Bifurcation curves for RM-model (2) and MB-model (54). Note that for RMmodel the expression are independent of ε. The arrow indicates the transition of the steady states that occurs when the parameter crosses the bifurcation point.

Table 4 :

 4 List of symbols used for RM-model with dimension with their meaning.the biomass of both populations have the same dimension.

	Symbol Meaning
	T	Dimensional time
	t	Dimensionless time
	X j	Dimensional state variable, indicated by j
	x j	Dimensionless state variable, indicated by j
	R	Intra-specific growth rate
	A	Attack rate
	F (X 1 ) Functional response, non-dimensional
	C	Efficiency, conversion yield, non-dimensional
	M	Mortality rate per unit of time
	k	Handling time

An other mathematical method to analyse a canard explosion is the blow-up technique [START_REF] Dumortier | Canard cycles and center manifolds[END_REF][START_REF] Dumortier | Geometric singular perturbation theory beyond normal hyperbolicity[END_REF][START_REF] Krupa | Relaxation oscillation and canard explosion[END_REF][START_REF] Kuehn | Multiple Time Scale Dynamics[END_REF]. This technique can be used to study the unfolding of the degenerated Hopf bifurcation where lim ε → 0. This will be the subject of a forthcoming paper.

In [START_REF] Muratori | Low-and high-frequency oscillations in three-dimensional food chain systems[END_REF][START_REF] Rinaldi | Slow fast limit-cycles in predator prey models[END_REF][START_REF] Rinaldi | Destabilizing factors in slow-fast systems[END_REF][START_REF] Hsu | Relaxation oscillation profile of limit cycle in perdator-prey system[END_REF][START_REF] Hek | Geometric singular perturbation theory in biological practice[END_REF][START_REF] Lobry | Effect of population size in a predator-prey model[END_REF] the RM-model predator-prey slow-fast model showed complex relaxation oscillation dynamics, however, no canard explosion was observed.

In [START_REF] Kooi | Aggregation methods in food chains with nutrient recycling[END_REF][START_REF] Auger | Aggregation methods in dynamical systems and applications in population and community dynamics[END_REF][START_REF] Kooi | Modelling the dynamics of traits involved in fighting-predators-prey system[END_REF] food chain systems were already studied where the small parameter, measuring the timescale disparity between the rate of changes in the prey (fast) and predator (slow), in the MB-model was introduced. Then the fast subsystem converges to a critical manifold of stable equilibria which yields an algebraic relationship between the state variable leading to a reduced (aggregated) lower dimensional (dae) system. Predator prey systems with fast oscillating migrations were studied in [START_REF] Poggiale | Fast oscillating migrations in a predator-prey model[END_REF] and with slow migrations in [START_REF] Marvá | Reduction of slow-fast periodic systems with applications to population dynamics models[END_REF] wherein reduction methods were proposed. In those papers, relaxation oscillations were not discussed.

In the MB-model, conservation of mass is obeyed and therefore it is more realistic than the RM-model where nutrients are not modelled explicitly. An open chemostat environment is assumed with inflowing nutrients and the outflow of all abiotic and biotic components.

Under specific situations the logistic growth rate of the prey is replaced by a prey growth model whereby nutrients are not only used for its own growth but also for its predator. This is a bottom-up effect in addition to the top-down effect of the prey consumption by the predator. The consumption occurs with a constant efficiency like in the alternative fast-slow version of the RM-model introduced in [START_REF] Hek | Geometric singular perturbation theory in biological practice[END_REF]Example 2.2]. The introduction of this constant efficiency in the fast-slow RM-model, however, leads to unrealistic equilibrium population sizes when ε becomes small. In the MB-model (54) these sizes remain realistic.

A constant efficiency, instead of the less realistic variable efficiency, was earlier used in [START_REF] Kooi | Aggregation methods in food chains with nutrient recycling[END_REF] in the fast-slow version of the RM-model and described by ( 2 The dynamics of the MB-model was analysed using a classical phase-space and bifurcation analysis approach where only equilibria and limit cycles occur in the whole parameter space. Calculations showed that in this more realistically and mechanistically underpinned model the complex fast-slow canard explosion does not occur.

A Derivation of the dimensionless RM-model

From [START_REF] Van Voorn | Combining bifurcation and sensitivity analysis for ecological models[END_REF] we recall (after some adjustments of the notation), that the classical RM-model is given as

where F (X 1 ) = X 1 (1 + AkX 1 ) -1 is the well-known Holling type II functional response, with X j the state variables, k handling time, A the attack rate, C a conversion efficiency, M the predator removal rate (mortality, maintenance, and harvesting), and T is time. A list of symbols is given in Table 4.

This model can be rescaled by using t = T R, x 1 = x 1 /K, and x 2 = X 2 /K. Note that in the last transformation for the predator biomass differs from [START_REF] Van Voorn | Combining bifurcation and sensitivity analysis for ecological models[END_REF] to let the efficiency C not disappear in the dimensionless formulation and to be able to deal with the time-scale difference in the fast-slow system the subject here.

For the special case b 1 = 4 at the Hopf bifurcation the results are shown in Fig. 19b.

The equilibrium E 2 of the full system (2) at point T in ( 14) is now for the reduced system restricted to the critical manifold not an equilibrium of the reduced system and there is no convergence to that point. Starting for x 2 (0) = 0.225 there is always convergence to the trivial zero solution E 0 which is here a global stable equilibrium of the reduced system.

In order to study the dynamics at the critical manifold (the parabola) further we plot [START_REF] Hsu | On global stability of a predator-prey system[END_REF] of the RM-model. In the top-subpanels the rate dx 1 /dt versus x 1 (t) while in the bottom-subpanel the dynamics x 2 (t) versus x 1 (t) on the critical manifold (see also Fig 19) Hopf bifurcation value b 1 = 4 the rate dx 1 /dt given in ( 21) is always negative and hence there is convergence always to the global stable zero solution. Note that the numerator in ( 21) is zero and this does not hold when b 1 = 4. However, for b 1 = 8 where the equilibrium is unstable the situation differs. Now when starting on the left side of point T the rate dx 1 /dt is positive and there is convergence toward the limit point T given in [START_REF] Ginoux | Flow curvature mathod applied to canard explosion[END_REF]. Similarly starting on the right-side the rate dx 1 /dt is negative and there is also convergence toward T . Now, the plot shows a vertical asymptote. Thus the rate is discontinuous at point T and this is due to the fact that dq/dx 1 (20) is zero at T . Calculations showed that the equilibrium point of the reduced system is reached in finite time. This is proved in the Appendix B.

In summary, the computational results show that point T is in the case of the reduced system not a simple tangent bifurcation. When the parameter equals the Hopf bifurcation value it is not even an equilibrium. Otherwise it is a limit point reached in finite time.

To support this the dynamics of the RM-model (2) on the parabola, the M 1 0 part of the f -nullcline (22b) is studied. From equations ( 20) and ( 21), one has

Since we have :

it follows that equation ( 57) is equivalent to

thus

where x 1 (T ) is the coordinate of the top of the parabola:

If we consider that the equilibrium point (

) is at point T , that is at the Hopf bifurcation point when b 1 = 4, it follows that: b

We can then express the time t T needed for starting from x 1 (0) to reach the equilibrium

But in this case point T is not an equilibrium. If one assumes that the equilibrium is on the right of point T , that is b 1 < 4, see Fig. 3a, (60) shows that the equilibrium is not reach in a finite time.

Finally, if we assume that the equilibrium x * 1 is on the left of point T that is b 1 > 4 and if we consider an initial condition between x * 1 and 1, see Fig.