CUBO, A Mathematical Journal Vol. 23, no. 03, pp. 343–355, December 2021 DOI: 10.4067/S0719-06462021000300343 Non-algebraic limit cycles in Holling type III zooplankton-phytoplankton models Homero G. D́ıaz-Maŕın1 Osvaldo Osuna2 1Facultad de Ciencias F́ısico-Matemáticas, Universidad Michoacana, Edif. Alfa, Ciudad Universitaria, C.P. 58040, Morelia, Michoacán, México. homero.diaz@umich.mx 2Instituto de F́ısica y Matemáticas, Universidad Michoacana, Edif. C-3, Ciudad Universitaria, C.P. 58040, Morelia, Michoacán, México. osvaldo.osuna@umich.mx ABSTRACT We prove that for certain polynomial differential equations in the plane arising from predator-prey type III models with generalized rational functional response, any algebraic solu- tion should be a rational function. As a consequence, limit cycles, which are unique for these dynamical systems, are necessarily trascendental ovals. We exemplify these findings by showing a numerical simulation within a system arising from zooplankton-phytoplankton dynamics. RESUMEN Probamos que para ciertas ecuaciones diferenciales poli- nomiales en el plano que aparecen a partir de modelos predador-presa de tipo III con respuesta funcional racional generalizada, toda solución algebraica debe ser una función racional. Como consecuencia, los ciclos ĺımite, que son únicos para estos sistemas dinámicos, son necesariamente óvalos trascendentes. Ejemplificamos estos resultados mostrando una simulación numérica para un sistema que aparece en la dinámica de zooplancton-fitoplancton. Keywords and Phrases: Predator-prey models, functional-response, Puiseux series, Newton polygon, limit cycles, invariant algebraic curve. 2020 AMS Mathematics Subject Classification: 34C25, 34M25, 34M35, 37N25, 92D25. Accepted: 07 July, 2021 Received: 08 January, 2021 ©2021 H. G. D́ıaz-Maŕın et al. This open access article is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License. http://cubo.ufro.cl/ http://dx.doi.org/10.4067/S0719-06462021000300343 https://orcid.org/0000-0002-2453-9049 344 H. G. Dı́az-Maŕın & O. Osuna CUBO 23, 3 (2021) 1 Introduction We consider predator-prey model u̇ = ru ( 1 − u K ) − vp(u), v̇ = v (−D + γp(u)) , (1.1) with functional response proposed in [13] p(u) = mun/(a + un). where the parameters are: the maximal feeding rate m; an affinity constant a, related to handling times, capture efficiencies, etc.; and the number of encounters n ≥ 1 a predator must have with a prey item before becoming maximally efficient at utilizing the prey item as a resource. According to [13], this last parameter is derived from an analogy with Michaelis-Menten equation for enzymatic kinetics. Here n measures the amount of ‘learning’ exhibited by the predator. For n > 1, this functional response has Holling type III, while for n = 1 it has Holling type II, that is why this functional is also called generalized functional response. Increasing the attack exponent 1 < n < 2 introduces the stability of simple consumer-resource population models, theoretical findings reveal that this increases biodiversity, see [14] and references therein. By fitting parameters, it is shown that n ≥ 2 appear in certain models in ecology, where predator free-space is a component of the habitat structure, see [1]. Other theoretical models of biological relevance consider the specific attack exponent n = 2, see [15, 18]. For 1 < n < 2 existence and uniqueness of limit cycles for predator-prey system (1.1) is proved in [16]. Existence and uniqueness for 0 ≤ n ≤ 1, n ≥ 2 also holds true under certain conditions, see [17]. Along this work we consider only integer values n ≥ 2. Existence of non-algebraic limit cycles for the Lotka-Volterra model were first exhibited by [12]. Since then, existence of trascendental ovals as limit cycles in system generalizing Lotka-Voltera models have been proved, see for instance [9, 5, 6]. Motivated by these results we explore this question for generalized functional responses. Our main result is contained in Theorem 2.1 which asserts that limit cycles can not be algebraic ovals in the case of Holing type III predator-prey models. The proof uses Puiseux series at infinity in the variable x. We estimate the number of branches of solutions given by the Puiseux series. We perform calculation of the leading term and prove that there exists at most one determination or branch of such series. To see this we show how each coefficient cn is completely specified by the parameters of the system. Thus, we conclude that any invariant algebraic curve must have at most degree one in y. Thus any algebraic invariant curve y(x), should be a rational function. For related works which also apply formal and Puiseux series to planar polynomial systems see [4, 7, 8]. CUBO 23, 3 (2021) Non-algebraic limit cycles in Holling type III ... 345 2 Rational functions as invariant algebraic curves Under a suitable change of variables, x = u, y = −v/m, and time reparametrization, ds dt = 1/(a + u2), system (1.1) becomes ẋ = x ( (r − 1/K)a + (r − 1/K)xn + xn−1y ) , ẏ = y (−D + (γm − D)xn) . (2.1) Thus we study the algebraic system: ẋ = x ( a0 + anx n + xn−1y ) , ẏ = y (b0 + bnx n) , an ̸= bn. (2.2) Notice that the axes x = 0, y = 0 are algebraic solutions of (2.1). Take the ODE defined by system (2.1) in the complex domain dy dx = y (b0 + bnx n) x (a0 + anxn + xn−1y) . (2.3) Solutions are Riemann surfaces immersed in Cx × Cy, where Cx ≃ C and poles of solutions correspond to values y = ∞ in the compactification Cy = Cy ∪ {∞} ≃ CP1. If we ask for the existence of algebraic solutions F(x, y) = 0 for F ∈ C[x, y], of the dynamical system (2.1). Then, such algebraic curve should be rational. Theorem 2.1. Suppose that the following conditions hold, a0 ̸= b0, an ̸= bn. (2.4) If there exists an invariant algebraic curve F(x, y) = 0 of equation (2.3) with x, y ∤ F(x, y), then degy F = 1. Therefore, any algebraic (possibly multivalued) solution should also be a rational (univalued) solution, y = ϕ(x), provided we exclude the trivial solution, y(x) ≡ 0. The following claim becomes of interest. Corollary 2.2. There can not exist algebraic limit cycles of the dynamical system (2.1) as a real vector field in R2, whenever conditions (2.4) hold true. For the proof of Theorem 2.1 we consider the Newton-Puiseux algorithm to describe explicitly the nature of solutions at the infinites x = ∞ and y = ∞. For further explanation of the Newton- Puiseux method for ODE, see [2, 10, 11]. The crucial step of the proof is to apply the following result. Theorem 2.3 (Theorem 1.4 in [3]). Let G(z, w) = 0 be an invariant algebraic curve, ∂wG ̸= 0 of the polynomial ODE P(z, w) dw dz − Q(z, w) = 0. (2.5) 346 H. G. Dı́az-Maŕın & O. Osuna CUBO 23, 3 (2021) Then degw G is at most the number of Puiseux series w(z) = c0z µ0 + ∞∑ l=1 clz l m0 +µ0, (2.6) solving (2.5), whenever the number of these series is finite. Here µ0 = l0/m0 with m0, l0 relatively prime integers m0 ≥ 0. Proof of Theorem 2.1. We proceed analyzing poles and algebraic branch points according to Painle- vé methodology, see [10, 11]. Notice that under the blow-up change of coordinates ξ = 1 x , equation (2.1) yields an equation at x = ∞ corresponding to ξ = 0 dy dξ = − y(a0ξ n + an) ξ(b0ξn + bn + ξy) , (2.7) At infinity the trivial solution y ≡ 0 yields a trivial solution which tends to ξ = 0. To find an expansion of non-trivial solutions along ξ = 0, with ξ = 1/x, in equation (2.7), we adopt the following Puiseux series expansion: y(ξ) = c0ξ µ0 + ∞∑ l=1 clξ l m0 +µ0, (2.8) where µ0 = l0/m0 and −1/µ0 is one of many possible slopes of the corresponding Newton polygon, and l0, m0 are relatively prime integers. For equation (2.7) the Newton polygon is a right-angled triangle whose only oblique side is the hypothenuse, see Fig. 1. y l 6 - . . .PP PP PP PPP� � � 2 1 2 3 4 5 1 n ξ l Figure 1: Newton polygon associated to the ODE (2.7) and used to calculate µ0. Circled vertices correspond to monomials appearing in By′ within the expression A(ξ, y) + B(ξ, y)dy dξ = 0. Therefore, the only slope to consider is −1/µ0 = 1. Accordingly, µ0 = −1 and c0(bn −an)−c20 = 0, with two possible roots: c0 = 0, bn −an ∈ C. If we make a direct substitution c0 = 0 of the Laurent expansion, ∑∞ l=0 clξ −1+l. This yields the trivial solution, y ≡ 0. We claim that the remaining value, c0 = bn − an ̸= 0 (2.9) CUBO 23, 3 (2021) Non-algebraic limit cycles in Holling type III ... 347 gives rise to a unique Laurent series of a simple pole at x = ∞, and therefore to just one branch of the Puiseux series. Indeed, under substitution ξ1 = ξ, y = c0ξ −1 1 + y1, we obtain a Newton polygon for (a0ξ n 1 + bn + ξ1y1)ξ1 dy1 dξ1 + (b0ξ n 1 + an)y1 + (an − bn)(a0 − b0)ξ n−1 1 = 0 (2.10) which has two possible slopes and corresponding values µ1 = 1, −1/(n − 1). See Fig. 2. n−1 l l 6 - . . . ..... PPPPPPPPP� � � PP PP PP PPP� � � 2 1 2 3 4 5 1 n ξ y l Figure 2: Newton polygon associated to the ODE (2.10) and used to calculate µ1. According to the algorithm given in [2], for positive slope −1/µ1, we choose as principal side, µ1 = n − 1. We have cn = (an − bn)(a0 − b0) (1 − n)bn − an . (2.11) In the following step, we have a principal side with µ2 = 2n − 1. See the corresponding Newton polygon used to calculate µ2 in Fig. 3. This determines c2n. Therefore, there is a unique determination for the Puiseux-Laurent series: y = c0ξ −1 + ∞∑ k=1 cknξ kn−1 + . . . . (2.12) By direct substitution of the Puiseux-Laurent series in eq. (2.7) we can also verify that the middle coefficients vanish, i.e. for each k = 0, 1, 2, . . . , we have cl = 0, ∀l = kn + 1, . . . , (k + 1)n − 1. Theorem 2.3 implies that degy F ≤ 1. Thus, under the hypothesis of Theorem 2.1 we conclude that y = ϕ(x) is a rational function which cannot contain an algebraic limit cycle because of the uniqueness of its determination with respect to x. 348 H. G. Dı́az-Maŕın & O. Osuna CUBO 23, 3 (2021) 1 l l - 6 . . . . � � aaaaaaaaaahhh hhh hhh hh1 y 2 3 4 n 2n − 1 l Figure 3: Newton polygon used to calculate µ2. Remark 2.4. We have chosen Puiseux-Laurent series of the form y = y(x) because there is a recognizable pattern in the successive Newton polygons, namely triangles with a moving low vertex. This yields a unique side with a unique slope. Therefore a unique µk yields a unique linear relation that allows us to compute all the coefficients ck. 3 On the degree with respect to x We may ask whether the degree in x for an invariant curve can be estimated with the same methods. Notice that expression (2.12) suggests that degx F = nk for some k ∈ N. We illustrate the difficulties to calculate an upper bound for degx F using the same techniques by considering n = 3. If we take, x = x(y) at y = ∞, then we may take the coordinate change y = 1 η . Thus system (2.1) becomes dx dη = a0xη + anx n+1η − xn−1 η2(b0 + bnxnη2) The corresponding Newton polygon is shown in Fig. 4. There are three posible cases for Puiseux-Laurent series x(η) = c0η µ0 + ∞∑ l=1 clη l m0 +µ0, corresponding to slopes −1/µ0 equal to 1, ∞, −(n − 1) which yield µ0 equal to −1, 0, 1n−1. No infinite values c0 ∈ R arise in each case. This can be verified as follows: (1) Case µ0 = −1. Under substitution η1 = η, x = c0η µ0 1 + x1, CUBO 23, 3 (2021) Non-algebraic limit cycles in Holling type III ... 349 x i - 6 B B B B B B B BB � � � n n+1 1 1 η i Figure 4: Newton polygon for x(y) at infinity used to calculate µ0 = 1 n−1. we regard the least degree coefficient. There exists a unique Puiseux series where c0 is determined by a linear relation cn0 (c0(bn − an) − 1) = 0, therefore there exists just one branch. (2) Case µ0 = 0. Corresponds to the trivial solution x ≡ 0 with c0 = 0. (3) Case µ0 = 1 n−1. Puiseux-Laurent series arise as c0 solve a relation: a0c0 + b0c0 n − 1 + cn0 = 0, therefore there exists n − 1 posible branches. Each branch corresponds to a (n − 1)−th root c0 = ( −a0 − b0 n − 1 )1/n−1 . If we choose µ0 = 1 n−1, then we get the ODE, A(η1, x1) + B(η1, x1) dx1 dη1 , with extended expression, A(3,0)η 3 1 +A(1,1)η1x1 + A(5/2,1)η 5/2 1 x1 + A(1/2,2)η 1/2 1 x 2 1 + A(2,2)η 2 1x 2 1 +A(0,3)x 3 1 + A(3/2,1)η 3/2 1 x1 + A(1,4)η1x 4 1 + dx1 dη1 × [B(1,1)η21 + B(5/2,1)η 7/2 +B(2,2)η 3 1x1 + B(3/2,3)η 5/2 1 x 2 1 + B(1,4)η 2 1x 4 1] = 0. 350 H. G. Dı́az-Maŕın & O. Osuna CUBO 23, 3 (2021) whose Newton polygon is shown in Fig 5. Circled vertices correspond to monomials appearing in Bx′1 within the ODE. Under the same assumptions of Theorem 2.1, if there exists an invariant algebraic curve F(x, y) = 0 of (2.3) with x, y ∤ F(x, y), then degx F has upper bound at least n − 1, provided we exclude the trivial solution, y(x) ≡ 0. This would require a0 + b0 n − 1 ̸= 0. (3.1) We still can not conclude that degx F ≤ n − 1, since the proof of this fact would require a suitable description of successive Newton polygons, as well as an effective calculation of the number of branches of the corresponding Puiseux-Laurent series. Two main difficulties arise: On one hand these Newton polygons may follow a complex pattern. On the other hand, we may have several different relations defining general coefficients ck, k > 0 requiring enough conditions so that there is a finite number of branches rather than a continuum. 4 l l ll 6 . . . . . . . . - A A A A H H HHA A A A A A AA� � x 1 2 3 η 1 2 3 l Figure 5: Newton polygon that determines µ1 = 2 with n = 3. In the second step we have the possibility to choose either µ1 = 2 or µ1 = 1/2. If we choose µ1 = 2. Then, the corresponding relations arising from the least degree terms in the substitution η1 = η2, x1 = c1η µ1 2 + x2 become, 4c1(4a0 − b0) = 8a3a20 + 4a 2 0b3 + 8a3a0b0 + 4a0b3b0 + 2a3b 2 0 + b3b 2 0, Therefore, we would require the additional condition 4a0 ̸= b0 (3.2) in order to be able to calculate c1 and thus have a finite number of branches. CUBO 23, 3 (2021) Non-algebraic limit cycles in Holling type III ... 351 According to [2, Lemma 2], in order to have a finite number of branches, i.e. a good side of the Newton Polygon in its terminology, it is sufficient that the following conditions hold for the high and low vertices of such a side, (a, b), (a′, b′), respectively: (1) B(a,b) ̸= 0 and A(a,b) B(a,b) /∈ Q≥µ1 = {q ∈ Q : q ≥ µ}, (2) A(a′,b′) + µB(a′,b′) ̸= 0, where A(a,b), B(a,b) refer to the coefficients for the monomials in the equation A + B dx1 dη1 = 0 associated to the vertex (a, b). In our concrete example, in Fig. 5 we have chosen µ1 = 2 because it corresponds to the slope −1/µ1 of the unique good side which has vertices (a, b) = (1, 1) and (a′, b′) = (3, 0). Calculations yield A(1,1) = 1 8 (16a0 + 12b0), B(1,1) = −b0, A(3,0) = −c40 ( a0 + b0 2 ) , B(3,0) = 0. Recall that under our conventions, B(3,0) = 0 is implied by the fact that the vertex (3, 0) is not circled. Conditions for a good side which are sufficient to have a finite number of branches read as follows: (1) 16a0+12b0 8b0 /∈ Q≥µ1 = {q ∈ Q : q ≥ µ1}. That is, either a0 b0 < 1 4 or a0 b0 ≥ 1 4 but a0 b0 /∈ Q. (3.3) (2) a0 + b0 2 ̸= 0. We recover condition (3.1). Notice that condition (3.3) is stronger than (3.2). In the following step we choose µ2 = 7/2 by considering the slopes of the Newton polygon shown in Fig. 6 with the unique good side which has vertex (a′, b′) on the η1−axis. Remark the increasing complexity of the polygon. Thus in the step k ≥ 2, we can always choose the good side largest negative slope −1/µk with vertex (a′, b′) = (a′, 0) on the ηk−axis and vertex (a, b) = (1, 1). But even if in each step k ≥ 2 we achieve linear relations to determine coefficients ck, we still can not conclude that there is a finite number of determinations. An additional calculation needs to be done, namely to verify that no other side in the Newton polygon, yield a continuous indetermination ck ∈ C. Those additional sides are not good. To illustrate this difficulty suppose that we do not choose the unique good side in Fig. 5. Suppose that on the contrary we choose the side with vertices (a, b) = (0, 3) and (a′, b′) = (1, 1) which is not good. Further calculation yields (a0 + 2b0)c0 = 0. 352 H. G. Dı́az-Maŕın & O. Osuna CUBO 23, 3 (2021) Therefore, either c0 = 0 or an indetermination c0 ∈ C arises whenever the following relation does or does not hold: a0 + 2b0 ̸= 0. (3.4) x l ll l l ll l l l. . . . . - . . . . . . . . . . . . 6 @ @XXXXXXX H H H H H H H H H H H H H HH� � 1 2 3 1 2 3 4 η4 5 6 7 8 9 l Figure 6: Newton polygon to determine µ2 = 7/2. Notice its increasing complexity with respect to Figs. 5 and 4. Summarizing, if we follow the same strategy, to find an effective upper bound for degx F requires further calculations and a detailed and complete description of the conditions that allow a finite number of branches. We leave it for a future work. 4 Zooplankton-phytoplankton dynamics We consider the dependence of a predator’s (zooplankton) grazing rate on prey (phytoplankton) is taken as that of Holling type III response as in [15], instead of type II as in (2.67), (2.68) in [18]. Suppose that phytoplankton grows in logistic form whereas the zooplankton predation by fish is neglected. We then get the following system u̇ = u (1 − u) − vu2 h + u2 , v̇ = γvu2 h + u2 − δv, (4.1) which yields a system similar to (2.1): ẋ = x (xy) , ẏ = y ( −δ + (γ − δ)x2 ) . (4.2) CUBO 23, 3 (2021) Non-algebraic limit cycles in Holling type III ... 353 The criterion exposed in [17] for existence and uniqueness of a limit cycle adapted to system (1.1) with m = 1 and n ≥ 2 states that (nD − (n − 2)γ) · n √ aD γ − D < (pD − (p − 1)γ)K, which for (4.1) yields 2δ √ hδ γ − δ < 2δ − γ. (4.3) For the specific choice of parameters: δ = 0.25, γ = 0.35, h = 0.01, condition (4.3) holds true. Therefore there exists a unique limit cycle. Indeed, numerical evidence for the existence of a limit cycle is given in Fig. 7. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.1 0.2 0.3 0.4 0.5 Figure 7: Solutions of system (4.1) with initial conditions (0.1, 0.05) and (0.2, 0.2) converge to a limit cycle. On the other hand, due to Corollary 2.2 this limit cycle can not be an algebraic curve. Indeed, for an algebraic invariant curve F(x, y) = 0, there is only one branch of a simple pole at x = ∞. The corresponding Puiseux-Laurent series of this branch is y = ∞∑ k=0 c2kξ 2k−1 = c0ξ −1 + c2ξ + c4ξ 3 + . . . , where ξ = 1/x. A straightforward calculation of the coefficients yields c0 = γ − δ = 0.1, c2 = δ = 0.25, c4 = 0 = c2k, k ≥ 2. Therefore, y = c0ξ −1 + c2ξ. Hence for an invariant algebraic curve, degy F = 1, and y = ϕ(x), should be rational with F(x, ϕ(x)) = 0. In this example degx F = 2. 354 H. G. Dı́az-Maŕın & O. Osuna CUBO 23, 3 (2021) References [1] D. Barrios-O’Neill, J. T. A. Dick, M. C. Emmerson, A. Ricciardi and H. J. MacIsaac, “Predator-free space, functional responses and biological invasions”, Functional Ecology, vol. 29, no. 3, pp. 377–384, 2015. [2] J. Cano, “An extension of the Newton-Puiseux polygon construction to give solutions of Pfaffian forms”, Ann. Inst. Fourier (Grenoble), vol. 43, no. 1, pp. 125–142, 1993. [3] M. V. Demina, “Novel algebraic aspects of Liouvillian integrability for two-dimensional poly- nomial dynamical systems”. Phys. Lett. A, vol. 382, no. 20, pp. 353–1360, 2018. [4] M. V. Demina, “Invariant algebraic curves for liénard dynamical systems revisited”, Appl. Math. Lett., vol. 84, pp. 42–48, 2018. [5] A. Ferragut and A. Gasull. “Non-algebraic oscillations for predator-prey models”, Publ. Mat., vol. 58, suppl., pp. 195–207, 2014. [6] J. Giné and M. Grau, “Coexistence of algebraic and non-algebraic limit cycles, explicitly given, using Riccati equations”, Nonlinearity, vol. 19, no. 8, pp. 1939–1950, 2006. [7] J. Giné and J. Llibre, “Strongly formal Weierstrass non-integrability for polynomial differential systems in C2”, Electron. J. Qual. Theory Differ. Equ., no. 1, pp. 1–16, 2020. [8] J. Giné and J. Llibre, “Formal Weierstrass nonintegrability criterion for some classes of poly- nomial differential systems in C2”, Internat. J. Bifur. Chaos Appl. Sci. Engrg., vol. 30, no. 4, 7 pages, 2020. [9] M. Hayashi, “On polynomial Liénard systems which have invariant algebraic curves”, Funkcial. Ekvac., vol. 39, no. 3, pp. 403–408, 1996. [10] E. Hille, Ordinary Differential Equations in the Complex Domain, Dover Publications, Inc., Mineola, NY, 1976. [11] E. L. Ince, Ordinary Differential Equations, Dover Publications, New York, 1944. [12] K. Odani, “The limit cycle of the van der Pol equation is not algebraic”, J. Differential Equations, vol. 115, no. 1, pp. 146–152, 1995. [13] L. A. Real, “The kinetics of functional response”, The American Naturalist, vol. 111, no. 978, pp. 289–300, 1977. [14] B. Rosenbaum and B. C. Rall, “Fitting functional responses: Direct parameter estimation by simulating differential equations”, Methods in Ecology and Evolution, vol. 9, no. 10, pp. 2076–2090, 2018. CUBO 23, 3 (2021) Non-algebraic limit cycles in Holling type III ... 355 [15] V. A. Ryabchenko, M. J. R. Fasham, B. A. Kagan and E. E. Popova, “What causes short-term oscillations in ecosystem models of the ocean mixed layer?”, Journal of Marine Systems, vol. 13, no. 1, pp. 33–50, 1997. [16] J. Sugie, “Uniqueness of limit cycles in a predator-prey system with Holling-type functional response”, Quart. Appl. Math., vol. 58, no. 3, pp. 577–590, 2000. [17] J. Sugie, R. Kohno, and R. Miyazaki, “On a predator-prey system of Holling type”, Proc. Amer. Math. Soc., vol. 125, no. 7, pp. 2041–2050, 1997. [18] R. K. Upadhyay and S. R. K. Iyengar, Introduction to Mathematical Modeling and Chaotic Dynamics, CRC Press, 2013. Introduction Rational functions as invariant algebraic curves On the degree with respect to x Zooplankton-phytoplankton dynamics