Modelling, singular perturbation and bifurcation analyses of bitrophic food chains B.W. Kooi, J.C. Poggiale To cite this version: B.W. Kooi, J.C. Poggiale. Modelling, singular perturbation and bifurcation analyses of bitrophic food chains. Mathematical Biosciences, Elsevier, 2018, 301, pp.93-110. ฀10.1016/j.mbs.2018.04.006฀. ฀halshs-01778648฀ HAL Id: halshs-01778648 https://halshs.archives-ouvertes.fr/halshs-01778648 Submitted on 23 Feb 2019 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. 1 2 Modelling, singular perturbation and bifurcation analyses of bitrophic food chains B.W. Kooi∗ and J.C. Poggiale∗∗, 3 ∗ ∗∗ Faculty of Science, VU Amsterdam, De Boelelaan 1085, 1081 HV Amsterdam, The Netherlands. Aix-Marseille University, UMR 7294 MIO OCEANOMED, Campus de Luminy, 163 Avenue de Luminy case 901, 13009 Marseille, France. 4 March 15, 2018 5 Abstract 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 Two predator-prey model formulations are studied: for the classical RosenzweigMacArthur (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 28 29 30 31 32 33 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. 36 Keywords: Aggregation methods; Asymptotic series expansion; Canard; Geometrical Singular Perturbation theory; Predator-prey models 37 Mathematics Subject Classification: 34E17 - 37G15 - 92D25 - 92D40 34 35 2 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 1 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 [44] 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 [45, 24]. 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 [19] it was proved that the coexistence equilibrium is globally stable and in [4] that the periodic solution is unique thus a globally stable limit cycle and in [31] 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” [43]. 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 [35, 42, 41, 5, 11]. 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 [41]). 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 [40]). The worm is the water nymph Nais elinguis, a oligochaete species. In a previous paper [25] 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 3 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 of the equilibria of the fast system but better approximations are obtained by asymptotic expansions [18, 16]. When the perturbation parameter becomes small the resulting limit cycle is called the relaxation oscillation, see [15], 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 [47, 10, 6, 3, 8, 28, 29, 2]. More recently, in [37, 16] geometric singular perturbation techniques [12] 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 [3, 2] 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. 2 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: a1 x2 dx1 = x1 (1 − x1 − ), dt 1 + b1 x1 c1 a1 x1 dx2 − d1 ) , = x2 ( dt 1 + b1 x1 104 105 106 107 108 109 110 111 112 113 (1a) (1b) with xi (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 a1 x1 /(1+b1x1 ) 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 a1 > 0 is the searching rate for the prey and b1 > 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 c1 > 0. The last term of (1b) is death rate d1 > 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. 4 114 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 c1 and death rate d1 both proportional to the positive perturbation parameter ε. For numerical studies we take that parameter a1 co-varies with b1 via a1 = 5/3 b1 , hence handling time is 3/5 and the values for parameters c = d = 1 for the RM-model and D = e1 = 1 for the MB-model. This is without loss of generality. parameter t τ x0 x1 x2 xr a1 b1 c1 = εc d1 = εd e1 D1 = εD ε ref.values [0, ∞) [0, ∞) [0, ∞) [0, ∞) [0, ∞) (0, ∞) a1 = 5/3 b1 3,4 or 8 c=1 d=1 1 D=1 [0, 1] Interpretation Fast time variable Slow time variable Nutrient density Prey biomass density Predator biomass density Nutrient concentration in reservoir Searching rate Handling time × searching rate Conversion efficiency Death rate Conversion efficiency Dilution rate Perturbation parameter 115 116 117 The set of equations analysed extensively in the literature that form a model with slow-fast dynamics reads a1 x2 dx1 = f (x1 , x2 , ε) = x1 (1 − x1 − ), dt 1 + b1 x1 a1 x1 dx2 − 1) , = εg(x1 , x2 , ε) = εx2 ( dt 1 + b1 x1 118 119 120 121 122 123 124 125 126 127 (2a) (2b) with xi ∈ R, i = 1, 2. We introduced c1 = εc and d1 = ε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 : R3 → R and g : R3 → 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 [42, 16]. 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 5 141 a small parameter. For, the efficiency c1 = ε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 (2). 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. 142 2.1 128 129 130 131 132 133 134 135 136 137 138 139 140 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 Existence and stability analysis of equilibria and limit cycles In model (2) there are only three free parameters, namely a1 , b1 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 b1 where ε = 1 and a1 co-varies with b1 via a1 = 5/3 b1 , is shown in Fig. 1 for parameter b1 as free parameter. The three relevant equilibria zero- E0 , boundary- E1 and interior equilibrium E2 and the limit cycle L2 are summarized in Table 2 and the bifurcation curves the transcritical bifurcation T C and Hopf bifurcation H in Table 3. In [19] it was proved that the coexistence equilibrium E2 is globally stable and in [4] that the periodic solution is unique thus a globally unique stable limit cycle L2 . 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 b1 . This phenomenon is related to the “paradox of enrichment” [43], because extinction due to stochastic fluctuations is likely. For the stability analysis of equilibria we need the Jacobian matrix evaluated at point (x1 , x2 ) ([27]), which reads J= 159  a1 x2 b1 a1 x1 a1 x2 + x1 ( (1+b − 1+b 1 − x1 − 1+b 2 − 1) 1 x1 1 x1 ) 1 x1 a1 x1 b1 a1 x1 εx2 ( 1+ba11 x1 − (1+b ) ε( − 1) 2 1+b1 x1 1 x1 )  . (3) In E0 the Jacobian matrix is   1 0 . J0 = 0 −ε 6 (4) TC H x2 1.0 0.5 0.0 x1 1.0 0.5 0.0 0 2 4 6 8 10 b1 Figure 1: One-parameter bifurcation diagram for b1 and a1 = 5/3 b1 of the RM-system (1) with ε = 1. The equilibria E1 (below T C) and E2 (between T C and H) as well as the maximum and minimum values for the limit cycle L2 (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 ε. Table 2: Equilibria of RM-model (2) and MB-model (54). Equilibria RM-model E0 = (0, 0) E1 = (1, 0) 1 , a1 −b1 −12 ) E2 = (x∗1 , x∗2 ) = ( a1 −b 1 (a1 −b1 ) MB-model E0 = (0, 0) E1 = (1, 0) a1 −b1 −1 1 E2 = (x∗1 , x∗2 ) = ( a1 −b , ) 2 1 ε(a1 −b1 ) +a1 −b1 7 System composition Extinction Prey-only Prey-predator Extinction Prey-only Prey-predator Table 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. Bifs. RM-model TC H MB-model TC H 160 161 a1 (b1 , ε) b1 (a1 = 5/3 b1 ) Interpretation a1 = b1 + 1 +1) a1 = b1b(b11−1 1.5 4 E1 → E2 E2 → L2 a1 = b1 + 1 √ 2εb2 +1+ 4b21 ε(ε+1)+1 a1 = 1 2ε(b1 −1) 1.5 b1 = √ 4ε+ 16ε2 +15ε (2ε) E1 → E2 E2 → L2 The eigenvalues are the diagonal elements and the E0 is always unstable. At the boundary equilibrium E1 the Jacobian matrix reads  −b1 −1  −1 − a11+b 1 . J1 = 1 −1 ) 0 ε( a1 −b a1 162 163 164 165 (5) The eigenvalues are again the diagonal elements and E1 is stable when a1 − b1 − 1 < 0 and unstable when a1 − b1 − 1 < 0. The equality gives the transcritical bifurcation T C where b1 = a1 − 1 or for the reference value a1 = 5/3 b1 we have b1 = 3/2. At the equilibrium E2 the Jacobian matrix reduces to J2 =  a1 (b1 −1)−b1 (1+b1 ) a1 ((a1 −b1 ) 1 −1 ε( a1 −b ) a1 −1 0  . (6) The complex pair of eigenvalues λ1,2 of the Jacobian matrix evaluated at the equilibrium E2 read λ1,2 (b1 , ε) = µ(b1 ) ± iω(b1 , ε) . 166 167 (7) We calculated the following real µ ∈ R and imaginary ω ∈ R parts as functions of parameter ε and b1 where a1 = 5/3 b1 −144 + 72b1 − 9b21 − 60b1 ε + 40εb21 , 25b21 3(b1 − 4) 1 − sgn ∆(b1 , ε)  ± |∆(b1 , ε)| , µ(b1 , ε) = 10b1 2 1 + sgn ∆(b1 , ε)  |∆(b1 , ε)| . ω(b1 , ε) = ± 2 ∆(b1 , ε) = (Tr J2 )2 − 4 det J2 = 8 (8) (9) (10) TC H 1 equilibrium ε 0.75 limit cycle ∆ < 0 spiral, focus 0.5 0.25 ∆ > 0 node 0 0 1 2 ∆ > 0 node 3 4 5 b1 6 7 8 9 10 Figure 2: Focus bifurcation where ∆ = 0 curve for ε vs b1 . 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 In Fig. 2 the curve where ∆(b1 , ε) = 0 is plotted in the (b1 , ε) 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 b1 = 3 or b1 = 8 the equilibria are foci while for ε = 0.01 the equilibria are nodes, stable and unstable, respectively. 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 b1 = 4 we have for ε > 0 a zero real part µ(4) = 0 and the√positive imaginary part equal to the determinant of the Jacobian matrix J2 : √ ω(4) = det J2 = 1/2 ε, for the two conjugated complex eigenvalues. The square root of the ratio µ(b1 )/ω(b1 ) measures the amplitude of the limit cycle that emerge for the Hopf bifurcation for b1 > 4. Evaluation of this ratio just above b1 = 4 gives that the oscillatory dynamics for values is as with the Hopf bifurcation till the canard explosion occurs. The derivative dµ/db1 (b1 )|b1 =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.3)] is satisfied. Using a similar procedure in Maple, [32] described in [31], the first Lyapunov coefficient ℓ1 is evaluated as ℓ1 = − 188 189 16 14ε + 1 . 75 ε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. 9 194 Furthermore, in the limiting case√ ε = 0 at the Hopf bifurcation where b1 = 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 [8, 9] which is beyond the scope of this paper. 195 2.2 190 191 192 193 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 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 b1 = 3, in the region between the transcritical T C and Hopf H bifurcation in Fig. 1 where equilibrium E2 is stable, • Fig. 4: for b1 = 4 at the Hopf bifurcation, • Right panels of Fig. 3: for b1 = 8 above the Hopf bifurcation in Fig. 1 where equilibrium E2 is unstable. In the top panels of Fig. 3 and Fig. 4 the f-nullcline where f (x1 , x2 , ε) = 0: x2 = (1 − x1 )(1 + b1 x1 )/a1 and the g-nullcline where g(x1 , x2 , ε) = 0: x1 = 1/(a1 − b1 ), are shown. Both are independent of ε. The graph of the f-nullcline is part of a parabola where it intersects the horizontal axis x2 = 0 at x1 = 1 and of the g-nullcline is just a vertical line through the equilibrium point E2 . 2.2.1 Stable interior equilibrium With b1 = 3, see Fig. 3a with ε = 1 and Fig. 3c with ε = 0.01, the equilibrium E2 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 d1 = 1: dx2 = −x2 , (12) dτ where differentiation is with respect to the slow time variable τ = εt. The solution crosses the critical point x2 = 1/a1 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 [42]. Eventually the system converges to the stable equilibrium E2 given in Table 2. 2.2.2 Hopf interior equilibrium The equilibrium E2 , see Table 2, at the Hopf bifurcation H coincides with the top of the f-nullcline parabola 10 1 1 0.8 0.8 ✷ ✷ x2 0.6 x2 0.6 0.4 0.4 ◦ 0.2 0.2 0 a ◦ 0 0 0.2 0.4 1 x1 0.6 0.8 1 b 0.8 0 0.2 0.4 1 x1 0.6 1 0.8 1 0.8 ✷ ✷ x2 0.6 x2 0.6 0.4 0.4 ◦ 0.2 0.2 0 c 0.8 ◦ 0 0 0.2 0.4 x1 0.6 0.8 1 d 0 0.2 0.4 x1 0.6 Figure 3: Phase-space analysis for system (2) describing the RM-model with b1 = 3 (left panels) and b1 = 8 (right panels) while a1 = 5/3 b1 . Top panels (a,b): ε = 1 and bottom panels (c,d): ε = 0.01. (a): Spiral stable equilibrium E2 , (b): stable limit cycle L2 , (c): node stable equilibrium E2 and (d): stable limit cycle L2 . 11 1 0.8 ✷ x2 0.6 0.4 ◦ 0.2 0 a 0 0.2 0.4 x1 1 0.6 0.8 1 1 0.8 0.8 ✷ ✷ x2 0.6 x2 0.6 0.4 0.4 ◦ 0.2 0 b ◦ 0.2 0 0 0.2 0.4 x1 0.6 0.8 1 c 0 0.2 0.4 x1 0.6 0.8 1 Figure 4: Phase-space analysis for system (2) describing the RM-model, b1 = 4 and a1 = 5/3 b1 , that is at the Hopf bifurcation point for three different values of ε: ε = 1 (a), ε = 0.1 (b) and ε = 0.01 (c). (x1 , x2 )T = 223  b1 − 1 (b1 + 1)2  . , 2b1 4a1 b1 (13) At the Hopf point we have the equilibrium E2 (x∗1 , x∗2 )H = (x1 , x2 )T =  3 15  . , 8 64 (14) 226 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 (14) at the top of the f-nullcline denoted by (13). 227 2.2.3 224 225 228 229 230 231 Unstable interior equilibrium With b1 = 8, see Fig. 3b,d equilibrium E2 is unstable and there is no convergence to E2 but to a stable limit cycle L2 . 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 L2 consists of four 12 0.3 g(x1 , x2 , 0) T ⋄ x2 M00 f (x1 , x2 , 0) M10 ◦ TC 0 0 0.2 0.4 0.6 0.8 1 x1 Figure 5: Fast (double arrow) and slow (single arrow arrow) dynamics for system (2) describing the RM-model with b1 = 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. 234 concatenated episodes. Two parts of the trajectory where the predator population changes slowly and two parts of the trajectory where the predator population is almost constant and the prey population changes fast (the almost horizontal parts of the solution orbit). 235 2.2.4 232 233 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 Degenerate phase curves With b1 > 4 the limits of the limit cycles L2 , the periodic relaxation oscillations, for lim ε → 0 are called degenerated phase curves shown in Fig. 5, see [14] 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 [42, 11] 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 (13). 13 252 253 254 255 256 257 258 259 260 261 262 Note that since b1 > 4 the equilibrium E2 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 b1 = 8 in both cases. However for the Hopf bifurcation b1 = 4 case the resulting dynamics for small ε ≪ 1 does not look like that of Fig. 5 and more importantly also not for b1 slightly above the Hopf bifurcation. Therefore we focus in the next sections on the dynamics for parameter values where the unstable equilibrium E2 is just above the Hopf bifurcation value. 3 Singular perturbation problem 266 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 [18, 21, 22] for introductions into perturbation analysis. 267 3.1 263 264 265 268 269 270 271 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: a1 x2 (0) dx1 = f (x1 , x2 (0), 0) = x1 (1 − x1 − ), dt 1 + b1 x1 dx2 =0. dt 272 273 274 (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: dx1 = f (x1 , x2 , ε) , dτ dx2 = εg(x1 , x2 , ε) . ε dτ ε 275 (15a) After substitution of ε = 0 we get 14 (16a) (16b) 0 = f (x1 , x2 , 0) , (17a) dx2 a1 x1 = g(x1 , x2 , 0) = x2 ( − 1) . dτ 1 + b1 x1 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 (17b) This differential algebraic equation (dae) is called the reduced system. The trajectory is either part of the vertical axis where x2 = 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 x1 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 (x1 , x2 , 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 (x1 , x2 , 0) = 0 is equivalent to x1 = p(x2 ) and we can substitute x1 by p(x2 ) in equation (17b). The result is the slow or reduced system dx2 (18) = g(p(x2 ), x2 ) . dτ The differentiation is with respect to τ and the smooth functions f and g are defined in (2). Then the slow system becomes 0 = f (p(x2 ), x2 , 0) , 292 293   (a1 − b1 )p(x2 ) − 1 dx2 , = x2 dτ b1 p(x2 ) + 1 294 295 296 x1 = p(x2 ) = An alternative method is to use instead of x1 = p(x2 ) of which the dynamics is described by the solution of (19), the inverse function x2 = q(x1 ) also derived from (17a). This gives the relationship x2 = q(x1 ) = 297   1  b1 − 1 ± (b1 + 1)2 − 4a1 b1 x2 . 2b1 (19) 1 (1 − x1 )(1 + b1 x1 ) , a1 and using (17b) the differential equation dq dx1 dx2 = g(x1 , q(x1 )) = , dτ dx1 dτ 15 (20) 298 and we get formally q(x1 )(a1 x1 − (1 + b1 x1 )) dx1 , = dτ (1 + b1 x1 ) dq/dx1  dq 1  = b1 (1 − 2x1 ) − 1 . dx1 a1 (21) 303 Note that this expression is zero at point T given in (13). Hence the denominator of model (21) 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 (21). Only when in the special case b1 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. 304 3.2 299 300 301 302 305 306 307 308 309 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 {(x1 , x2 )|f (x1 , x2 , 0) = 0, x1 ≥ 0, x2 ≥ 0} consist of two types of critical manifolds,see Fig. 5 M00 = (x1 , x2 )|x1 = 0, x2 ≥ 0 1 M01 = (x1 , x2 )|x2 = (1 − x1 )(1 + b1 x1 ), x1 , x2 ≥ 0 . a1 310 311 312 313 314 316 317 318 319 (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 = F (X, Y, ε) , dt dY = εG(X, Y, ε) , dt dε =0, dt 315 (22a) (23a) (23b) (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 DF set D, all the eigenvalues of 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); 16 320 321 322 323 324 325 326 327 328 329 330 331 ii) the graph of p is invariant under the flow associated to the original differential system (23); 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 M00 and M10 are normally hyperbolic and there exists ε0 such that for 0 < ε < ε0 , there are locally invariant manifolds M0ε and M1ε except in the neighborhood of point T (x1 , x2 ) = (x1 , x2 ) and in the neighborhood of the intersection between M00 and M10 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 M1ε can be approximated by asymptotic expansions in ε. It can be described as a graph {(x1 , x2 )|x2 = q(x1 , ε), x1 ≥ 0, x2 ≥ 0} . 332 This manifold is invariant when the following equality holds dx2 dx1 dq(x1 , ε) dx1 dx2 = = , dt dx1 dt dx1 dt 333 (26) Then (23) gives with x2 = q(x1 , ε) the invariance condition  x1 (a1 − b1 ) − 1  ∂q(x1 , ε)  a1 q(x1 , ε)  x1 1 − x1 − = εq(x1 , ε) , ∂x1 1 + b1 x1 1 + b1 x1 335 (25) yields with equation (23) and x2 = q(x1 , ε):  x1 (a1 − b1 ) − 1  ∂q(x1 , ε) dx1 = εq(x1 , ε) . ∂x1 dt 1 + b1 x1 334 (24) (27) or using 1 + b1 x1 > 0    ∂q(x1 , ε)  x1 (1 − x1 )(1 + b1 x1 ) − a1 q(x1 , ε) = εq(x1 , ε) x1 (a1 − b1 ) − 1 . ∂x1 336 3.2.1 337 The following asymptotic expansion in ε is introduced: (28) Asymptotic expansion q(x1 , ε) = q0 (x1 ) + εq1 (x1 ) + ε2 q2 (x1 ) + . . . , 17 (29) 338 hence dq2 ∂q dq0 dq1 = +ε + ε2 + ... . ∂x1 dx1 dx1 dx1 (30) Substitution into (28) gives 339   dq2 dq1 dq0 +ε + ε2 + . . .)x1 (1 − x1 )(1 + b1 x1 ) − a1 (q0 (x1 ) + εq1 (x1 ) + ε2 q2 (x1 ) + . . .) dx1 dx1 dx1    (31) = ε(q0 (x1 ) + εq1 (x1 ) + . . .) x1 (a1 − b1 ) − 1 . ( 340 Gathering orders of ε results for O(1) and assuming x1 > 0 in: q0 (x1 ) = 341 342 343 (1 − x1 )(1 + b1 x1 ) , a1 344 346 gives   ∂q(x1 , ε) x1 a1 (q0 − q(x1 , ε)) = εq(x1 , ε) x1 (a1 − b1 ) − 1 , ∂x1 348 349 350 351 (x1 (a1 − b1 ) − 1) . dq0 −a1 x1 dx 1 (33) (34) At b1 = 4 the numerator and denominator are both zero. we have x1 = (b1 − 1)/(2b1 ) and dq0 /dx1 = 0 but also since it is a equilibrium x1 (a1 − b1 ) − 1 = 0. For O(ε2 ) in (33) gives q2 (x1 ) =q1 (x1 ) 347 (32) At b1 = 4 we have x1 = (b1 − 1)/(2b1 ) and hence dq0 /dx1 = 0. For O(ε) and using an updated form of (28) q1 (x1 ) = q0 (x1 ) 345 dq0 b1 − 1 − 2x1 b1 = . dx1 a1 dq1 xa dx1 1 1 + x1 (a1 − b1 ) − 1 dq0 −a1 x1 dx 1 . (35) At point T we have x1 = x1 = (b1 − 1)/(2b1 ) where dq0 /dx1 = 0. Consequently, at that point the denominator in the expression for qi (x1 ), i > 0 is zero. Therefore the coefficients qi (x1 ), i > 0 are unbounded when the numerator is not equal zero. Only when the parameter b1 = 4 is at the Hopf bifurcation the numerator is zero and the coefficients qi (x1 ) remain finite. 18 352 353 354 355 356 357 358 359 3.2.2 Asymptotic expansion in phase space The expression for q0 describes the critical manifold M10 . This expression is the inverse (when it exists) of p(x2 ) in (19). The voluminous expressions for the higher order qi , i > 1 coefficients obtained by equating the O(εi ) terms on the left- and right-hand side of the invariance condition (28), are not given here but are available using Maple, [32]. This yields the approximation of the perturbed slow manifold M1ε For ε = 0 the limit (16b) prescribes the singular slow flow on M10 with x2 = q0 (x1 ) given by (32) dx1 a1 q0 (x1 ) = x1 (1 − x1 − ). dt 1 + b1 x1 360 361 362 For sufficiently small non-zero ε ≪ 1 the flow on the perturbed slow manifold M1ε can be approximated by inserting x2 = q(x1 , ε) with q(x1 , ε) given by (29). In order to simulate the model we solve  dx1 a1 q(x1 , ε)  , = x1 1 − x1 − dt 1 + b1 x1 363 364 365 (36) (37) with properly chosen the initial values. The results are shown in Fig. 6 for b1 = 3, where the equilibrium E2 is stable. They show that in this case the solution of the original model on M1ε is already well approximated by the second order approximation when ε = 0.1. 0.3 ◦ 0.25 x2 0.2 0.15 ✷ 0.1 0.05 0 0 0.2 0.4 x1 0.6 0.8 1 Figure 6: Results for the original system (2) (solid) describing the RM-model (29) with b1 = 3 where a1 = 5/3 b1 and initial conditions x1 = 0.86038 and x2 (0) = 0.102646. 366 367 368 369 In Fig. 7a and Fig. 7b the graph of the function q(x1 , ε) is shown for b1 = 3 and b1 = 8, respectively. These results show that the asymptotic expansion for ε > 0 is only locally a good approximation for M1ε but fails at the top of the parabola point T . 19 0.3 0.3 ◦ 0.25 0.25 0.2 x2 x2 0.2 0.15 0.1 0.1 0.05 0.05 0 a ◦ 0.15 0 0 0.2 0.4 x1 0.6 0.8 1 0 b 0.2 0.4 x1 0.6 0.8 1 0.3 0.25 ◦ x2 0.2 0.15 0.1 0.05 0 c 0 0.2 0.4 x1 0.6 0.8 1 Figure 7: Second-order approximation of M1ε results with (a): b1 = 3 (b): b1 = 8 and (c): b1 = 4 where a1 = 5/3 b1 for the RM-model (29). Solid: with ε = 0.01. Dotted with ε = 0.1. 370 371 372 373 374 At the Hopf bifurcation In Fig. 7c the graph of the function q(x1 , ε) is shown for b1 = 4. While the solution converges to the neutral stable equilibrium E2 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 x1 = x∗1 = x1 : lim q(x1 , ε) = ∞ and x1 ↓x1 lim q(x1 , ε) = −∞ . x1 ↑x1 (38) 376 In the next section we extend the asymptotic expansion in ε by varying parameter b1 in addition to ε in order to repair this unwanted property. 377 3.3 375 378 379 380 381 382 Canard explosion In this section we analyse the Canard dynamics that occurs for b1 values just above the Hopf bifurcation at b1 = 4 similar to the analysis performed in [2]. Other papers on canards are [10, 6, 8, 11, 2]. 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. 20 383 384 385 386 387 We start with an exploration of this oscillatory behaviour for b1 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 b1 = 4. The results are shown for ε = 0.01 in Fig. 8. H x2 1.0 0.5 0.0 1.0 x1 0.5 0.0 3.99 4 4.01 4.02 b1 4.03 4.04 4.05 Figure 8: One-parameter bifurcation diagram for b1 and a1 = 5/3 b1 of the RM-system (1) with ε = 0.01. 388 389 390 391 392 393 394 395 Unexpectedly the amplitude of the limit cycle increases sharply when b1 passes a value just above b1 = 4.04. In Fig. 9 parameter ε is varied continuously for a number values of b1 just above the Hopf bifurcation point. Using (8) these results support the analytical expression for the amplitude of the  unique stable limit cycles emerging from the Hopf bifurcation for small ε values being: µ(b1 )/ω(b1 ). The continuation of the curves calculated with AUTO [7] failed for small ε. This is due to the part of the cycle close to the vertical axis where x1 is small, see Fig. 3bd. In order to avoid this dynamics we study now an augmented system dx1 a1 x2 = δ + f (x1 , x2 , ε) = δ + x1 (1 − x1 − ), dt 1 + b1 x1 a1 x1 dx2 = εg(x1 , x2 , ε) = εx2 ( − 1) , dt 1 + b1 x1 396 397 398 399 400 (39a) (39b) where δ is a small allochthonous input rate of the prey population. Addition of this extra term removes the transcritical bifurcation at x2 = 1/a1 because it is structurally unstable with respect to such a perturbation. Fig. 10 is a similar diagram as Fig. 9 where δ = 0.0001 instead of δ = 0. Note that the Hopf bifurcation occurs at values slightly different form b1 = 4, namely at b1 = 4.000711364 and this is taken into account in what follows. The 21 x2 0.25 0.125 0.0 b1 = 4.01 4.02 x1 1.0 4.03 4.04 4.05 0.5 0.0 0 Figure 0.002 0.004 0.006 0.008 0.01 0.012 0.014 ε 9: One-parameter bifurcation diagram for ε with various b1 = 4.01, 4.02, 4.03, 4.04, 4.0401, 4.0407, 4.05 values where a1 = 5/3 b1 of the RM-system (1). For b1 > 4 the minimum and minimum populations values during the stable limit cycle ate show. Also the two dashed curves are shown for b1 = 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. 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 results in Fig. 10 imply that continuation is possible toward very low ε values. The results 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 b1 just below and above the sudden changes: b1 = 4.04019 and b1 = 4.04061 where ε = 0.01. These results show how the unique stable limit cycle changes shape very abruptly at b1 ≈ 4.0403 where b1 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 b1 = 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 22 x2 0.25 0.125 0.0 4.005 4.01 x1 1.0 4.02 4.03 4.04 0.5 0.0 0 Figure 0.002 0.004 0.006 0.008 0.01 0.012 0.014 ε 10: One-parameter bifurcation diagram for ε with various b1 = 4.00125, 4.0025, 4.005, 4.01, 4.02, 4.03, 4.04 values where a1 = 5/3 b1 of system (39) with prey input rate δ = 0.0001. The parameter continuation is here successful to very small values of ε. 421 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. 422 3.3.1 420 423 424 425 426 427 428 429 430 431 432 Asymptotic expansion and canard explosion Following [13] 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 x2 = q(x1 , ε) evaluated at q0 (x1 ) (see (32)) but also in the bifurcation parameter b1 = b1 (ε) and consequently a1 (ε) = 5/3 b1 (ε) evaluated at the Hopf bifurcation point b10 = 4 with a01 = 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(x1 , ε) and the expansion is evaluated at r(x1 = x1 , ε = 0). When b1 = b10 we have r(x1 , ε) = q(x1 , ε). Similar to the invariance condition of the perturbed slow manifold M1ε we derive now the adapted equation that replaces (28) where we substitute a1 (ε) = 5/3 b1 (ε):    ∂r  x1 (1 − x1 )(1 + b1 (ε)x1 ) − 5/3 b1 (ε)r(x1 , ε) = εr(x1 , ε) 2/3 x1 b1 (ε) − 1 . ∂x1 433 (40) The following extended asymptotic expansion in ε for r(x1 , ε) is introduced as follows: x2 = r(x1 , ε) = r0 (x1 ) + εr1 (x1 ) + ε2 r2 (x1 ) + . . . , 23 (41) 0.3 0.25 x2 0.2 0.15 0.1 0.05 0 0 0.2 0.4 0.6 0.8 1 x1 Figure 11: Phase-space diagram with ε = 0.01. Solid lines: stable limit cycle for two values of b1 = 4.04019 (small cycle) and b1 = 4.04061 (big cycle) while a1 = 5/3 b1 of the RM-system (1). Dashed line: extended asymptotic expansion r(x1 , ε) of M1ε (41) where b1 = 4.0403. 0.3 0.237 0.25 ◦ ◦ 0.2 x1 x1 ε = 0.01 0.15 ε = 0.008 0.227 ε = 0.009 0.1 0.05 ε = 0.011 0 0.217 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 x2 x2 Figure 12: Phase-space diagram with ε = 0.008, 0.009, 0.01, 0.011. Solid lines: stable limit cycle for two values ε = 0.008, 0.009 (big cycle) and two values ε = 0.01, 0.011 (small cycle) of the RMsystem (1). Dashed line: extended asymptotic expansion r(x1 , ε) of M1ε (41) where b1 = 4.0403. 24 0.6 434 435 436 437 438 and dr2 ∂r dr0 dr1 = +ε + ε2 + ... , ∂x1 dx1 dx1 dx1 (42) b1 (ε) = b10 + εb11 + ε2 b12 + ε3 b13 + ε4 b14 . . . , (43) whereby b1 (ε) is given by where rj and b1j , 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(1) terms yields: r0 = 439 440 441 443 444 445 446 447 448 449 450 451 (44) With b10 = 4 this gives the zero order approximation of M1ε equal to M10 . Using (44) and substitution of (41), (42) and (43) in (40) and equating the resulting first order O(ε) terms, yields: r1 (x1 ) = 442 (1 − x1 )(1 + b10 x1 ) . 5/3 b10 3(1 − x1 )(1 − b10 + 2x1 b10 )b11 x1 b10 − x1 b210 − 3b10 + 2x21 b310 ) . 5b210 (1 + 2x1 b10 − b10 )x1 (45) It appears that the terms with b11 , namely 3x1 b11 (b10 − 1 − 2x1 b10 ) = 0, are already eliminated when they are evaluated at the Hopf bifurcation support point with b10 = 4 and a10 = 20/3 at the equilibrium E2 with x1 = x1 = x∗1 = 1/(a10 − b10 ). At that point both the numerator and the denominator are zero but limx1 ↓x1 r1 (x1 ) = limx1 ↑x11 r1 (x1 and this shows that function r1 (x1 ) is continuous for all b11 . From this point a recursive procedure can be followed. One by one, ri , i ≥ 2 is determined by taking the ith order term in the invariance condition (40) equal to zero and thereafter the term bth 1(i−1) by the condition that this term is continuous at x1 = x1 = (b10 − 1)/(2b10 ). 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 r2 = (1 − x1 ) (−288x41 + 108x31 )b12 + (72x41 − 27x31 )b211 + (96x31 + 84x21 )b11 + 256x31 − 128x21 − 112x1 − 16 . 960x31 (8x1 − 3) (46) 452 453 At the support point T being the Hopf bifurcation point with b10 = 4 and x1 = x1 given by (14) we have (−288x41 + 108x31 )b12 = 0. Furthermore the denominator is zero b10 = 4. The 25 454 455 456 457 458 459 460 461 462 expression for b11 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 x2 = q(x1 , ε) at q(x1 , ε) for coefficient q1 given in (34) where this expression was also unbounded. Now, because we have more freedom we can take b11 so that the expression stays bounded at this point. The expression b11 = 100/27 is obtained by substitution of the equilibrium value for x1 = x1 = 3/8 given by (14) into (45). Further recursion gives higher order approximations, again first rn with bn−1 then rn+1 with bn and so on. We calculated the following fourth order approximation b1 (ε) = b10 + εb11 + ε2 b12 + . . . , 100 58700 80536900 171270040300 b1 (ε) = 4 + ε + ε2 + ε3 + ε4 . 27 2187 177147 14348907 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 (47) In Fig. 11 besides the shape of the limit cycles for two values b1 = 4.04019 (small limit cycle) and b1 = 4.04061 (big limit cycle) also the result of the extended asymptotic expansion r(x1 , ε), where ε = 0.01, (41) with b1 = 4.0403 is shown. The approximation for M1ε 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 x1 changes sign and the rate becomes fast. This occurs for the smallest b1 = 4.04019-value. For the big limit cycle at b1 = 4.04061, from the separation point the trajectory along the cycle continuous to move toward the vertical axis. The asymptotic expansion for M1ε with b1 ≈ 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 b1 ≈ 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(x1 , ε) 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(x1 , ε) (41) is divergent as shown in Fig. 13 where the coefficients b1i as function of i is depicted. It was shown in [39, 38] 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 b1 parameter value is taken from Fig. 10. Fortunately, the result presented are in agreement. 26 4 2 log10 b1iεi 0 -2 -4 -6 -8 0 5 10 15 20 25 i 30 35 40 45 50 Figure 13: Coefficients b1i as function of i given by (47) with ε = 0.01. 0.05 b1 − 4.000711364 0.04 0.03 0.02 0.01 0 0 0.004 0.008 ε 0.012 Figure 14: Distance from parameter value b1 where canard explosion occurs from b1 = 4.000711364 where the Hopf bifurcation occurs with allochthonous prey input δ = 0.0001 as function of ε. Solid line is graph of the truncated expression (43) and the dots taken from Fig. 10. 27 487 488 489 490 491 492 4 The MB bitrophic food chain model This section presents the Monod chemostat model [34, 45, 24, 26]. 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 x0 (t) denote the density of the nutrient, and xi (t) ∈ R+ , t ≥ 0, i = 1, 2 the biomass densities of prey and predator, respectively. The scaled version of the Monod model reads dx0 = (xr − x0 )D1 − a0 x0 x1 , dt a1 x1 x2 dx1 , = a0 x0 x1 − D1 x1 − dt 1 + b1 x1 a1 x1 x2 dx2 − D1 x2 . = e1 dt 1 + b1 x1 493 494 495 496 497 498 499 500 501 t≥0. 505 (49) (50) This gives dH = −D1 H , dt   a1 x1 x2  dx1  , = H + xr − x1 − x2 x1 − D1 x1 + D1 dt 1 + b1 x1   a1 x1 dx2 = D1 x2 −1 . dt 1 + b1 x1 504 (48c) So, the asymptotic behavior of system (48) is bounded by xr . It is possible to decouple this system by the introduction of the function H(t) = x0 (t) + x1 (t) + x2 (t) − xr , 503 (48b) where xr is the concentration of nutrient in the reservoir and D1 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 e1 is the ratio of these two terms. For simplicity we assume e1 = 1 and a0 = 1. It can be shown that all solutions system (48) starting in the non-negative cone eventually lie in the set Ω = x = (x0 , x1 , x2 ) ∈ R3+ : x0 + x1 + x2 ≤ xr . 502 (48a) (51a) (51b) (51c) When furthermore H = 0 for the asymptotic dynamics we can study the two dimensional predator-prey system 28   a1 x2 dx1 = x1 xr − x1 − x2 − D1 − D1 , dt 1 + b1 x1  a1 x1  dx2 −1 . = D1 x2 dt 1 + b1 x1 506 507 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 (53) In order to be able to make a clear comparison with the RM-model formulations possible we use xr = 1 + D1 . This gives with D1 = ε   a1 x2 dx1 , = x1 1 − x1 − x2 − ε dt 1 + b1 x1  a1 x1  dx2 = εx2 −1 . dt 1 + b1 x1 508 (52) (54a) (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 x1 (1 − x1 ) is replaced by the expression x1 (1 − x1 − x2 ) with one extra term namely that of the predator biomass, x2 . In absence of the predator, x2 = 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 −D1 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 [45] (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 x2 term in the expression for the prey growth rate. 4.1 Existence and stability analysis of equilibria and limit cycles For the analysis of model is the chemostat environment we refer to [45]. 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 ε. Firstly we calculate by continuation of the parameter b1 the bifurcation diagram shown in Fig. 15. This diagram looks very much the same as Fig. 1 for the RM model. The main difference is that the Hopf bifurcation occurs at a somewhat higher b1 value. 29 TC H x2 1.0 0.5 0.0 x1 1.0 0.5 0.0 0 2 4 6 8 10 b1 Figure 15: One-parameter bifurcation diagram for b1 and a1 = 5/3 b1 of the MB-system (54) with ε = 1. 534 4.2 Phase-space analysis 540 Figure 16 displays the simulation results for the MB-model for three values of ε where b1 = 3. There is convergence to a stable equilibrium, similar as we found for the RMmodel in Fig. 3a,c. For b1 = 8 the results are shown in Fig. 17. These results differ much from the RMmodel in Fig. 3b,d case where the equilibrium E2 was unstable. Here this holds true for ε = 1 (Fig. 17a) but for smaller values the equilibrium becomes stable (Fig. 17b,c). 541 4.3 535 536 537 538 539 The degenerate phase point 549 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 x0 (0) + x1 (0) and the predator population remains constant x2 (0). Hence, the equilibrium E2 is neutral stable. The degenerate phase curve is just this point which is an equilibrium point x1 together with the initial predator size x2 (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. 550 4.4 542 543 544 545 546 547 548 551 552 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 b1 (whereby a1 = 5/3 b1 ), parameter ε is the second variable. 30 1 0.8 ✷ x2 0.6 0.4 0.2 ◦ 0 0 a 0.2 0.4 x1 1 0.6 0.8 1 1 0.8 0.8 ✷ ✷ x2 0.6 x2 0.6 ◦ 0.4 0.2 0.2 0 b ◦ 0.4 0 0 0.2 0.4 x1 0.6 0.8 1 c 0 0.2 0.4 x1 0.6 0.8 1 Figure 16: Phase-space analysis for system (54) describing the MB-model with a1 = 5/3 b1 , b1 = 3, for three different values of ε: ε = 1, ε = 0.1 and ε = 0.01. 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 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 (x2 = 0) equilibrium E1 is also the same. For the MBmodel the Hopf bifurcation H terminates at the origin where b1 → ∞ and ε → 0. There is a stable equilibrium E2 in almost the whole b1 > 4 range up to lim b1 ↑ ∞ while in the RM-model there is a stable limit cycle L2 . 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. 5 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 [17] and Michaelis-Menten kinetics. In this paper we compared two fast-slow versions of predator-prey models: the RM31 1 0.8 ✷ x2 0.6 0.4 0.2 ◦ 0 0 a 0.2 0.4 x1 1 0.6 0.8 1 0.8 0.8 ◦ ✷ 0.6 ✷ 0.6 x2 x2 ◦ 0.4 0.4 0.2 0.2 0 b 1 0 0 0.2 0.4 x1 0.6 0.8 1 c 0 0.2 0.4 x1 0.6 0.8 1 Figure 17: Phase-space analysis for system (54) describing the MB-model with a1 = 5/3 b1 , and b1 = 8, that is at the Hopf bifurcation point for three different values of ε: ε = 1, ε = 0.1 and ε = 0.01. 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 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 [35, 42, 41] the dynamics of this bitrophic systems has been studied intensively. In [11] this subject was also studied focusing on the dynamics close to the critical manifold M00 (Fig. 5) (part of the vertical axis where the prey is absent) and related to a delayed (transcritical) bifurcation phenomenon. Recently in [16] 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 M1ε slow manifold close to the parabolic critical manifold M10 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 [2] a similar method has been used: here the advantage is that all calculations are done in Maple, [32] with rational numbers. Direct application, however, showed that the delayed transcritical bifurcation dynamics close to the critical manifold M00 leads to an disturbing effect because the prey population 32 1/8 1/5 1/4 1/3 1 0.8 stable limit cycle stable equilibrium ε 0.6 0.4 0.2 HMB HRM TC 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 1/b1 Figure 18: Two-parameter bifurcation diagram with ε and 1/b1 as free parameters. The expressions are given in Table 3. The and MB-model the same. For lim ε ↓ 1 they differ essentially. whole b1 > 4 range up to lim b1 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 transcritical bifurcation curve T C is for both models RM-model ε = 1 the point of the Hopf bifurcations differ slightly but for In the MB-model there is a stable equilibrium E2 in almost the ↑ ∞ while there is a stable limit cycle L2 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 [46] 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 b1 and perturbation parameter ε values where a canard explosion occurs (see Fig. 11). 33 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 An other mathematical method to analyse a canard explosion is the blow-up technique [8, 9, 29, 30]. 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 [35, 42, 41, 20, 16, 11] the RM-model predator-prey slow-fast model showed complex relaxation oscillation dynamics, however, no canard explosion was observed. In [26, 1, 23] 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 [36] and with slow migrations in [33] 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 [16, 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 [26] in the fast-slow version of the RM-model and described by (2) [16, Example 2.1]. 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. References [1] P. Auger, R. Bravo de la Parra, J-C. Poggiale, E. Sanchez, and L. Sanz. Aggregation methods in dynamical systems and applications in population and community dynamics. Phys Life Rev, 5(2):79–105, 2008. [2] M. Brons. An iterative method for the canard explosion in general planar systems. Discrete and continuous dynamical systems, 250:77–83, 2013. [3] M. Canalis-Durand. Formal expansion of van der pol equation canard solutions are gevrey. In E. Benoı̂t, editor, Dynamic Bifurcation, pages 28–39. Springer, 1990. [4] K.-S. Cheng. Uniqueness of a limit cycle for predator-prey system. SIAM Journal on Mathematical Analysis, 12:541–548, 1981. 34 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 [5] B. Deng. Food chain chaos with canard explosion. Chaos, 14(4):1083–1092, 2004. [6] M. Diener. The canard unchained or how fast/slow dynamical problems bifurcate. The Mathematical Intelligencer, 6:38–49, 1984. [7] E. J. Doedel and B. Oldeman. Auto 07p: Continuation and bifurcation software for ordinary differential equations. Technical report, Concordia University, Montreal, Canada, 2009. [8] F. Dumortier and R. Roussarie. Canard cycles and center manifolds, volume 121 of Memoires of the American Mathematical Society. American Mathematical Society, AMS, Providence, RI, USA, 1996. [9] F. Dumortier and R. Roussarie. Geometric singular perturbation theory beyond normal hyperbolicity. In C. K. R. T. Jones and A. I. Khibnik, editors, Multiple Time Scale Dynamical Systems, volume 122 of IMA, pages 29–64. Springer-Verlag, Berlin, 2000. [10] W. Eckhaus. Relaxation oscillations including a standard chase on french ducks. In Asymptotic analysis II, volume 985 of Lecture Notes in Mathematics, pages 449–494, Berlin, 1983. Springer-Verlag. 666 [11] C. Lobry F. Campillo. Effect of population size in a predator-prey model. Ecol Model, 246:1–10, 2012. 667 [12] N. Fenichel. Geometric singular perturbation theory. JDE, 31:53–98, 1979. 665 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 [13] J-M. Ginoux and J. Llibre. Flow curvature mathod applied to canard explosion. J. Phys. A Math. Theor., 44(46):465203, 2016. [14] J. Guckenheimer. Normal Forms, Bifurcations and Finiteness Problems in Differential Equations, volume 137 of NATO Sci. Ser. II Math. Phys. Chem., chapter Bifurcations of relaxation oscillations, pages 295–316. Kluwer, Dordrecht, The Netherlands, 2004. [15] J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields, volume 42 of Applied Mathematical Sciences. SpringerVerlag, New York, 2 edition, 1985. [16] G. Hek. Geometric singular perturbation theory in biological practice. J Math Biol, 60:347–386, 2010. [17] C. S. Holling. Some characteristics of simple types of predation and parasitism. Canadian Entomologist, 91:385–398, 1959. [18] F. C. Hoppensteadt. Analysis and Simulation of Chaotic Systems. Applied Mathematical Sciences. Springer-Verlag, Berlin, 1993. [19] S.-B. Hsu. On global stability of a predator-prey system. Math Biosci, 174(1-2):1–10, 1978. [20] S.-B. Hsu and J. Shi. Relaxation oscillation profile of limit cycle in perdator-prey system. Discrete and continuous dynamical systems series B, 11(4):893–911, 2009. 35 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 [21] C. K. R. T. Jones. Geometric singular perturbation theory. Dynamical Systems, 1609:44–118, 1995. [22] J. Kevorkian and J.D. Cole. Multiple Scale and Singular Perturbation Methods, volume 114 of Applied Mathematical Sciences. Springer-Verlag, Berlin, 1995. [23] B. W. Kooi. Modelling the dynamics of traits involved in fighting-predators–prey system. J Math Biol, 71:1575–1605, 2016. [24] B. W. Kooi, M. P. Boer, and S. A. L. M. Kooijman. Consequences of population models on the dynamics of food chains. Math Biosci, 153(2):99–124, 1998. [25] B. W. Kooi, J. C. Poggiale, and P. Auger. Aggregation methods in food chains. Math. Comp. Mod., 27(4):109–120, 1998. [26] B. W. Kooi, J. C. Poggiale, P. Auger, and S. A. L. M. Kooijman. Aggregation methods in food chains with nutrient recycling. Ecol Model, 157(1):69–86, 2002. [27] M. Kot. Elements of Mathematical Ecology. Cambridge University Press, Cambridge, 2001. [28] M. Krupa and P. Szmolyan. Geometric analysis of the singularly perturbed fold. In J. K. R. T. Christopher and A Khibnik, editors, Multiple-Time-Scale Dynamical Systems, volume 122 of The IMA Volumes in Mathematics and its Applications, pages 89–116. Springer, 2001. [29] M. Krupa and P. Szmolyan. Relaxation oscillation and canard explosion. Journal of Differential Equations, 174:312–368, 2001. [30] C. Kuehn. Multiple Time Scale Dynamics, volume 191 of Applied Mathematical Sciences. Springer-Verlag, New York, 2015. 709 [31] Yu. A. Kuznetsov. Elements of Applied Bifurcation Theory, volume 112 of Applied Mathematical Sciences. Springer-Verlag, New York, 3 edition, 2004. 710 [32] Maple. Maple software. Maplesoft, Waterloo, Ontario, Canada, 2008. 708 711 712 713 714 715 716 717 718 719 720 721 722 [33] M. Marvá, J-C. Poggiale, and R. Bravo de la Parra. Reduction of slow-fast periodic systems with applications to population dynamics models. Mathematical Models and Methods in Applied Sciences, 22(10), 2012. [34] J. Monod. Recherches sur la croissance des cultures bactériennes. Hermann, Paris, 1942. [35] S. Muratori and S. Rinaldi. Low- and high-frequency oscillations in three-dimensional food chain systems. SIAM J Appl Math, 52:1688–1706, 1992. [36] J. C. Poggiale and P. Auger. Fast oscillating migrations in a predator-prey model. Mathematical Models & Methods in Applied Sciences (M3AS), 6(2):217–226, 1996. [37] J. C. Poggiale, P. Auger, F. Cordoleani, and T. Nguyen-Huu. Study of a virus-bacteria interaction model in a chemostat:application of geometrical singular perturbation theory. Philos T Roy Soc B, 367:4685–3428, 2009. 36 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 [38] J. P. Ramis. Les développements asymptotiques après poincaré : continuit et. . . divergences. La gazette des mathématiciens, 134:17–36, 2012. [39] J. P. Ramis. Poincaré et les développements asymptotiques. mathématiciens, 133:33–72, 2012. La gazette des [40] C. H. Ratsak, S. A. L. M. Kooijman, and B. W. Kooi. Modelling the growth of an oligocheate on activated sludge. Wat. Res., 27(5):739–747, 1993. [41] S. Rinaldi and A. Gragnani. Destabilizing factors in slow–fast systems. Ecol Model, 180:445–460, 2004. [42] S. Rinaldi and S. Muratori. Slow fast limit-cycles in predator prey models. Ecol Model, 61(3-4):287–308, 1992. [43] M. L. Rosenzweig. Paradox of enrichment: destabilization of exploitation ecosystems in ecological time. Science, 171:385–387, 1971. [44] M. L. Rosenzweig and R. H. MacArthur. Graphical representation and stability conditions of predator-prey interactions. Am Nat, 97:209–223, 1963. [45] H. L. Smith and P. Waltman. The Theory of the Chemostat. Cambridge University Press, Cambridge, 1994. 740 [46] J. Sugie and Y. Saito. Uniqueness of limit cycles in a Rosenzweig-MacArthur model with prey immigration. SIAM J Appl Math, 72(1):299–316, 2012. 741 [47] B. Van der Pol. On relaxation oscillations. Philosophical Magazine, 7:978–992, 1926. 739 742 743 744 [48] G. A. K. van Voorn and B. W. Kooi. Combining bifurcation and sensitivity analysis for ecological models. Eur. Phys. J. Special Topics, 226:2101–2118, 2017. A Derivation of the dimensionless RM-model From [48] we recall (after some adjustments of the notation), that the classical RM-model is given as dX1 = RX1 (1 − X1 /K) − AF (X1 )X2 , dT dX2 = (CAF (X1 ) − M) X2 , dT 745 746 747 748 749 750 751 752 (55a) (55b) where F (X1 ) = X1 (1 + AkX1 )−1 is the well-known Holling type II functional response, with Xj 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, x1 = x1 /K, and x2 = X2 /K. Note that in the last transformation for the predator biomass differs from [48] 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. 37 Table 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 t Xj xj R A F (X1 ) C M k Dimensional time Dimensionless time Dimensional state variable, indicated by j Dimensionless state variable, indicated by j Intra-specific growth rate Attack rate Functional response, non-dimensional Efficiency, conversion yield, non-dimensional Mortality rate per unit of time Handling time The non-dimensional model then reads 753 dx1 a1 x1 x2 , = x1 (1 − x1 ) − dt 1 + b1 x1 c1 a1 x1 dx2 − d1 ) , = x2 ( dt 1 + b1 x1 754 where a1 = AK/R, b1 = kAK, c1 = C and d1 = M/R. 755 B 756 757 758 759 760 761 762 763 764 765 766 767 768 (56a) (56b) The slow dynamics on the parabola part of the f nullcline 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 (12). 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 x1 (τ ) while (20) is used to get the associated predator population size x2 (τ ). Figure 19a gives the results for the stable equilibrium E2 , b1 = 3 where a1 = 5/3 b1 = 5 case with initial condition x2 (0) = 0.225. This value is below point T , where x2 = 0.25. In this case there are two valid initial points given by (19). Starting at the largest x1 prey value on the right-hand parabolic branch there is convergence toward the equilibrium E2 . But at the lowest prey value on the left-hand branch there is convergence to the trivial zero-solution where x1 = 0 and x2 = 0. 38 0.3 0.3 ✷•◦ 0.25 ✷ 0.25 ✷ ✷ a x2 0.15 x2 ✷ 0.15 0.1 0.1 0.05 0.05 0 ◦ 0 0.2 0.4 x1 0.6 0.8 1 b 0 ◦ 0 0.3 0.3 0.25 0.25 0.2 ✷◦ 0.15 0.2 0.2 • ✷ ◦✷ 0.15 0.1 0.1 0.05 0.05 0 c ◦ 0.2 x2 x2 0.2 0.4 x1 0.6 • 0.8 1 0.8 1 ✷ 0 0 0.2 0.4 x1 0.6 0.8 1 d 0 0.2 0.4 x1 0.6 Figure 19: Phase-space analysis for slow system (19) of the RM-model. (a): At the top-panel left b1 = 3 and a1 = 5/3 b1 , with initial values x2 = 0.225 and x2 = 0.256. The equilibrium E2 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 b1 = 4 and a1 = 5/3 b1 , with initial values x2 = 0.225. The equilibrium E2 is unstable. (c) and (d): The lower panel with b1 = 8 and a1 = 5/3 b1 the two initial values x2 = 0.15 and x2 = 0.16. The associated initial value of the x1 (0) is for x2 = 0.15 below the unstable equilibrium E2 value. For x2 = 0.16 the associated initial value of the x1 (0) is above the unstable equilibrium E2 value. 769 770 771 772 773 774 775 776 777 778 779 780 781 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 E2 . On the other hand for the left-hand branch it is in the basin of attraction of the zero solution. Figure 19c,d were calculated with parameter values b1 = 8 and a1 = 5/3 b1 = 40/3 where the equilibrium E2 of the original system is unstable. In Fig. 19c starting at the lowest prey value and x2 (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 E2 of the full model. Note that the vector field is not defined at the top T . In Fig. 19d the initial value is x2 (0) = 0.16 where both simulations terminate at a limit point T of the reduced system. Hence, the unstable equilibrium E2 of the full system is a separatrix between the two equilibrium points of the reduced system, being limit point T in (13) and the zero point E0 (x1 , x2 ) = (0, 0). 39 782 783 784 785 786 787 For the special case b1 = 4 at the Hopf bifurcation the results are shown in Fig. 19b. The equilibrium E2 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 x2 (0) = 0.225 there is always convergence to the trivial zero solution E0 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 dx1 /dt versus x1 (t) in Fig 20 in addition to x2 (t) versus x1 (t) in the phase space. For the 10.0 ◦ 0.0 -2.0 0.25 dx1/dt dx1/dt 2.0 ✷ ◦ -10.0 0.25 ✷ 0.125 • ◦✷ ✷ x2 x2 0.125 a • 0.0 0.0 ◦ 0 0.2 0.4 x1 0.6 0.8 1 b 0.0 ◦ 0 0.2 0.4 x1 0.6 0.8 1 Figure 20: Phase-space analysis for slow system (19) of the RM-model. In the top-subpanels the rate dx1 /dt versus x1 (t) while in the bottom-subpanel the dynamics x2 (t) versus x1 (t) on the critical manifold (see also Fig 19). (a): The left-panel b1 = 4 and a1 = 5/3 b1 , with initial value x2 = 0.225 where the equilibrium E2 is stable. (b): The right-panel with b1 = 8 and a1 = 5/3 b1 the initial value is x2 = 0.16 whereby the associated initial value of the x1 (0) is above the unstable equilibrium value. 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 Hopf bifurcation value b1 = 4 the rate dx1 /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 b1 = 4. However, for b1 = 8 where the equilibrium is unstable the situation differs. Now when starting on the left side of point T the rate dx1 /dt is positive and there is convergence toward the limit point T given in (13). Similarly starting on the right-side the rate dx1 /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/dx1 (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 M10 part of the f -nullcline (22b) is studied. From equations (20) and (21), one has (1 − x1 )((a1 − b1 )x1 − 1) dx1 = . dτ b1 − 2b1 x1 − 1 40 (57) 804 Since we have : b2 + b1 − a1 b1 + a1 b1 + 1 1 1 −2b1 x1 + b1 − 1 + 1 = , (1 − x1 )((a1 − b1 )x1 − 1) b1 + 1 − a1 1 − x1 b1 + 1 − a1 (a1 − b1 )x1 − 1 (58) 805 it follows that equation (57) is equivalent to b1 + 1 b21 + b1 − a1 b1 + a1 dx1 = (b1 + 1 − a1 )dτ , + 1 − x1 (a1 − b1 )x1 − 1 806 thus x1 (T ) x1 (0) 807 808 809 810 811 b1 + 1 b21 + b1 − a1 b1 + a1  + dx1 = 1 − x1 (a1 − b1 )x1 − 1 813 814 815 816 817 818 819 0 (b1 + 1 − a1 )dτ , (60) b1 − 1 . 2b1 1 ) is at point T , that is at the a1 − b1 Hopf bifurcation point when b1 = 4, it follows that: b21 + b1 − a1 b1 + a1 = 0. We can then express the time tT needed for starting from x1 (0) to reach the equilibrium E2 at point T If we consider that the equilibrium point (x∗1 = 1 (1 + b1 )/(2b1 )  (1 + b1 ) ln( ) . a1 − b1 − 1 1 − x1 (0) (61) 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 b1 < 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 b1 > 4 and if we consider an initial condition between x∗1 and 1, see Fig. 3d, point T attracts the trajectory. It is, however, not an equilibrium in the usual sense because equation (57) does not vanish, it is actually not well defined. Nevertheless this point is reached in a finite time according to (60), and the time needed to reach this point t∗T is t∗T 820 T where x1 (T ) is the coordinate of the top of the parabola: x1 (T ) = x1 = tT = 812 (59)  1 (1 + b1 )/(2b1 ) A/(2b1 ) A = ln( (1 + b1 ) ln( )+ ) , a1 − b1 − 1 1 − x1 (0) a1 − b1 (a1 − b1 )x1 (0) − 1 where A = a1 b1 − a1 − b21 − b1 > 0 under the above mentioned conditions. 41 (62)