This is a preprint of: “Phase portraits of the Selkov model in the Poincaré disc”, Jaume Llibre, Arefeh Nabavi, Discrete Contin. Dyn. Syst. Ser. B, vol. 27(12), 7607–7623, 2022. DOI: [10.3934/dcdsb.2022056] PHASE PORTRAITS OF THE SELKOV MODEL IN THE POINCARÉ DISC JAUME LLIBRE1 AND AREFEH NABAVI2 Abstract. In this paper we classify the phase portraits in the Poincaré disc of the Selkov model ẋ = −x + ay + x2 y, ẏ = b − ay − x2 y, in function of its parameters a, b ∈ R. We determine the regions in the parameter plane with biological meaning. 1. Introduction and statement of the main results In this paper we consider the Selkov model of glycolisis which is given by the following cubic differential system (1) ẋ = −x + ay + x2 y, ẏ = b − ay − x2 y, where x and y are the concentrations, and a and b are kinetic parame- ters. As usual the dot denotes derivative with respect to the time t. For more biological meanings about this system see the papers [10, 11], and the references quoted there. For other related models see for instance [3]. Our purpose is to classify the global phase portraits of system (1) in the Poincaré disc, in function of its parameters, and determine the regions of the parameters with biological meaning. Roughly speaking the Poincaré disc is the closed disc D2 centered at the origin of R2 of radius one, its interior is identified with R2 and its boundary, the circle S1 is identified with the infinity of R2 , because in R2 we can go to the infinity in many directions as points has the circle S1 . Then polynomial vector fields X, as the ones defined by system (1), can be extended to analytic vector fields p(X) defined in 2010 Mathematics Subject Classification. Primary 34C05. Secondary 34C07, 37C08. Key words and phrases. Selkov model, Phase portrait, Poincaré compactification. 1 2 J. LLIBRE AND A. NABAVI D2 . See subsection 2.2 for precise definitions related with the Poincaré compactification. The symmetry (x, y, b) → (−x, −y, −b) leaves invariant the differen- tial system (1), so we only need to investigate this system for the value of the parameter b > 0. Our main result is the following one. Theorem 1. The global phase portraits of the Selkov systems (1) in the Poincaré disc are topologically equivalent to one of the 9 phase portraits of Figure 1. A B C D E F G H I Figure 1. Phase portraits of system (1) with a ∈ R and b > 0 in the Poincaré disc. The shaded areas correspond to the initial conditions of the orbits having a finite final evolution, so these are the initial conditions with biological meaning. In the phase portraits E and I the final behaviour is a stable limit cycle, while in the phase portraits F , G and H the final behaviour is a singular point. SELKOV MODEL 3 a R1 b2 b -0.25 b1 δ2 δ1 R2 β Figure 2. The regions R1 and R2 where the differential system (1) has limit cycles computed numerically. The proof of Theorem 1 is given in section 3. The proof of Theo- rem 1 is done modulo the following conjecture about the areas in the parameter space {(b, a) : b ≥ 0, a ∈ R} where the limit cycles of the differential system (1) exist. Conjecture The differential system (1) has a unique stable limit cycle in the region R1 of Figure 2, and a unique unstable limit cycle in the region R2 of Figure 2. We shall prove that the limit cycle in the Selkov system appears from a Hopf bifurcation when √ the values of the parameters (b, a) are in the 2 curve a = (−1 − 2b + 1 + 8b2 )/2, and disappears in a graphic having a part at infinity, for this last claim we only have numerical evidences. In [2] it was done the following conjecture: The Selkov system (1) with b = 1 has a unique limit cycle when a ∈ (a1 , 0) where a1 ∈ (−0.036, −0, 037). Moreover, this limit cycle is stable and it is born in a graphic with a part at infinity when a = a1 and ends in a Hopf bifurcation at the singular point (1, 1/(1 + a)) when a = 0. Recently this conjecture has been proved by Chen and Tang in [4]. In the paper [2] there is the classification of the phase portraits in the Poincaré disc of the Selkov system (1) with b = 1 when a varies. In this particular subfamily of Selkov systems appear the phase portraits of Theorem 1 with the exception of the phase portraits A, G and I of Figure 1. Some definitions and basic results necessary for proving Theorem 1 are given in the section 2. The proof of Theorem 1 done in section 3 is divided in five subsections, in the first we study the local phase portraits of the finite singular points, in the second the local phase 4 J. LLIBRE AND A. NABAVI portraits of the infinite singular points, in the third the limit cycles, in the fourth the local phase portraits of all finite and infinite singular points in the Poincaré disc, and in the last the global phase portraits in the Poincaré disc. 2. Preliminary results 2.1. Singular points. Consider the following polynomial differential system in R2 (2) ẋ = P (x, y), ẏ = Q(x, y). The degree of the polynomial differential system (2) is the maximum of the degrees of the polynomials P (x, y) and Q(x, y). A point (x0 , y0 ) ∈ R2 is the singular point of system (2), when P (x0 , y0 ) = Q(x0 , y0) = 0. If all the eigenvalues of the Jacobian matrix   ∂ ∂ P (x0 , y0 ) P (x0 , y0 )  ∂x ∂y    M =   ∂ ∂  Q(x0 , y0 ) Q(x0 , y0 ) ∂x ∂y have nonzero real parts, the singular point (x0 , y0 ) is hyperbolic, and if one of the eigenvalues is zero the singular point is semi-hyperbolic. If both of the eigenvalues of the matrix M are zero and the matrix is not identically zero, the singular point (x0 , y0 ) is nilpotent. If M is identically zero, the singular point (x0 , y0 ) is linearly zero. For determining the local phase portraits of the hyperbolic, semi- hyperbolic and nilpotent singular points we can use Theorems 2.15, 2.19 and 3.5 of [5], respectively. For linearly zero singular points, we need to use a special change of variables called blow-ups, see [1] for details. 2.2. Poincaré compactification. In this section we provide a short summary about the Poincaré compactification. This method allows to study the behaviour of a planar polynomial differential system in a neighborhood of the infinity. See Chapter 5 of [5] for additional details on the Poincaré compactification. Let X = (P (x, y), Q(x, y)) be a polynomial vector field of degree d. Consider the Poincaré sphere S2 = {y = (y1 , y2 , y3 ) ∈ R3 : y12 +y22 +y32 = SELKOV MODEL 5 1} and its tangent plane at the point (0, 0, 1) which is identified with R2 , we assume that in this R2 we have defined our vector field X. Define the central projection f : R2 → S2 which sends every point x ∈ R2 , to the two intersection points f (x) of the straight line passing through the point x and (0, 0, 0) with the sphere S2 . It is clear that the equator S1 = {y ∈ S2 : y3 = 0} of the sphere S2 corresponds with the infinity of R2 . The differential Df sends the vector field X on R2 into a vector field X′ on S2 \ S1 . The vector field X′ essentially is formed by two symmetric copies with respect to the origin of the vector field X, one in the northern hemisphere of S2 and the other in the southern hemisphere. The vector field X′ can be extended analytically to a vector field on the whole sphere S2 multiplying X′ by y3d . We denote this new vector field on S2 by p(X), and it is called the Poincaré compactification of the polynomial vector field X on R2 . For studying the dynamics of X in the neighborhood of the infinity, we must study the dynamics of p(X) near S1 . The sphere S2 is a 2– dimensional manifold so we need to know the expressions of the vector field p(X) in the local charts (Ui , φi ) and (Vi , ψi ), where Ui = {y ∈ S2 : yi > 0}, Vi = {y ∈ S2 : yi < 0}, φi : Ui −→ R2 and ψi : Vi −→ R2 for i = 1, 2, 3, with φi (y) = −ψi (y) = (ym /yi , yn /yi ) for m < n and m, n 6= i. In the local chart (U1 , φ1 ) the expression of p(X) is        d 1 u 1 u d+1 1 u u̇ = v −uP , +Q , , v̇ = −v P , . v v v v v v In the local chat (U2 , φ2) the expression of p(X) is        d u 1 u 1 d+1 u 1 u̇ = v P , − uQ , , v̇ = −v Q , , v v v v v v and in the local chart (U3 , φ3 ) the expression of p(X) is u̇ = P (u, v), v̇ = Q(u, v). The expressions for p(X) in the local chart (Vi , ψi ) is the same than in the local chart (Ui , φi ) multiplied by (−1)d−1 for i = 1, 2, 3. The equator S1 is invariant by the vector field p(X) and all the singular points of p(X) which lie in this equator are called the infinite singular points of X. If y ∈ S1 is an infinite singular point, then −y is also an infinite singular point and these two points have the same (respectively opposite) stability if the degree of vector field is odd (respectively even). 6 J. LLIBRE AND A. NABAVI The image of the northern hemisphere of S2 onto the plane y3 = 0 under the orthogonal projection π(y1 , y2 , y3) = (y1 , y2 ) is called the Poincaré disc denoted by D2 . Since the orbits of p(X) on S2 are sym- metric with respect to the origin of R3 , we only need to consider the flow of p(X) in the closed northern hemisphere. So we project by π the phase portrait of p(X) on the northern hemisphere onto the Poincaré disc D2 , and we denote the compactified vector field in the Poincaré disc by π(p(X)). We shall present the phase portraits of the polynomial differential systems (1) in the Poincaré disc. The points at the infinity have the coordinates (u, 0) in the local charts Ui and Vi with i = 1, 2. Note that it is sufficient to study the infinite singular points on the local chart U1 , and then check the singularity of the origin of the local chart U2 . 2.3. Topological equivalence of two polynomial vector fields. The two polynomial vector fields X1 and X2 on R2 are topologically equivalent if there exists a homeomorphism on the Poincaré disc D2 which preserves the infinity S1 and sends the trajectories of the flow of π(p(X1 )) to the trajectories of the flow of π(p(X2 )), preserving or reversing the orientation of all the orbits. A separatrix of the Poincaré compactification π(p(X)) is an orbit at the infinity S1 , or a finite singular point, or a limit cycle, or an orbit on the boundary of a hyperbolic sector at a finite or an infinite singular point. The set of all separatrices of π(p(X)) is a closed (see [8]) and we denote it by ΣX . An open connected component of D2 \ ΣX is a canonical region of π(p(X)). The separatrix configuration of π(p(X)) is the union of an orbit of each canonical region with the set ΣX , and it is denoted by Σ′X . We denote by S (respectively R) the number of separatrices (re- spectively canonical region) of a vector field π(p(X)). We say that two separatrix configurations Σ′X1 and Σ′X2 are topolog- ically equivalent if there is a homeomorphism h : D2 −→ D2 such that h(Σ′X1 ) = Σ′X2 . The following theorem of Markus [7], Neumann [8] and Peixoto [9] allows to investigate only the separatrix configuration of a polynomial differential system in order to determine its global phase portrait. Theorem 2. The phase portraits in the Poincaré disc of two compact- ified polynomial vector fields π(p(X1 )) and π(p(X2 )) with finitely many SELKOV MODEL 7 separatrices are topologically equivalent if and only if their separatrix configurations Σ′X1 and Σ′X2 are topologically equivalent. 3. Proof of Theorem 1 3.1. Finite singular points. Case 1: Assume that a 6= −b2 . Then for any values a ∈ R and b > 0 the unique singular point of system (1) is p = (b, b/(a + b2 )). Linearizing this system at p we obtain the matrix  2  b −a 2  a + b2 a+b  J =  . −2b 2  2 −a − b a + b2 Then the determinant and trace of J are a + b2 and −(a − b2 + (a + b2 )2 )/(a + b2 ), respectively. For determining the local phase portrait of the singular point p we also need the eigenvalues of the matrix J, which are p (a − b2 + (a + b2 )2 ) ± (a(1 + a) + (2a − 1)b2 + b4 )2 − 4(a + b2 )3 − . 2(a + b2 ) In order to control the signs of (a(1+a)+(2a−1)b2 +b4 )2 −4(a+b2 )3 , of the determinant and of the trace of the matrix J we define the curves: 1 √ p √ α1 : a = (1 − 2 2b − 2b2 ∓ 1 − 4 2b), 2 1 √ p √ α2 : a = (1 + 2 2b − 2b2 + 1 + 4 2b), 2 1 √ p √ α3 : a = (1 + 2 2b − 2b2 − 1 + 4 2b), 2 1 √ β : a = (−1 − 2b2 + 1 + 8b2 ), 2 γ : a = −b2 . shown in Figure 3. It is easy to investigate the stability of the singular point p on these curves and in the regions (1)-(6) bounded by them, see Figure 3 and Table 1. For the values of the parameters (a, b) on the points √ of the curve β the trace of the matrix J is zero and its determinant is ( 8b2 + 1−1)/2 > 0. So the singular point p of the differential system (1) is a weak focus for these values of the parameter. Then the stability or instability of this weak focus is controlled by its Lyapunov constants, for more details 8 J. LLIBRE AND A. NABAVI a α2 α1 (1) (2) (3) b2 b -0.25 γ (4) (6) (5) β α3 Figure 3. Regions and curves separating the different local phase portraits of the finite singular point. Curve Region finite singular point α1 , α2 (1), (2) stable node β (b 6 b2 ) (3) stable focus β (b > b2 ) (4) unstable focus α3 (5) unstable node γ empty (6) saddle Table 1. We provide the local phase portrait of the finite singular point when the parameters of the differential system (1) are on the curves or in the regions of Figure 3. on these constants see Chapter 4 of [5]. Computing the first Lyapunov constants we get p√ √  √  π 8b2 + 1 − 1 4 −4b2 + 8b2 + 1 + 2 b2 + 8b2 + 1 + 1 V3 = − √ . 32 2b4 q √ This Lyapunov constant V3 = 0 only when b = b2 = 3/4 + 3/2, and if b ∈ (0, b2 ) then V3 < 0, while if b > b2 then V3 > 0. Therefore, if b < b2 the weak focus p is stable, and if b > b2 the weak focus p is unstable. When b = b2 then the next non–zero Lyapunov constant is π V5 = − q √  = −0.0380507785... 4 3 71 + 41 3 SELKOV MODEL 9 Hence for b = b2 we have a stable weak focus at p. Case 2: Assume that a = −b2 . Then the differential system (1) has no finite singular points. 3.2. Infinite singular points. Now we shall study the infinite singu- lar points of the polynomial vector field X. From subsection 2.2 we need to use the expressions of the compactified vector field p(X) in the local charts U1 and U2 . Thus the Poincaré compactification of system (1) in the local chart (U1 , φ1) is u̇ = −u − u2 + (1 − a)uv 2 + bv 3 − au2 v 2 , v̇ = −v(u − v 2 + auv 2 ). This system has two singular points on v = 0 which are (−1, 0) and (0, 0). The eigenvalues of the Jacobian matrix at (−1, 0) are 1, 1 and this point is a hyperbolic unstable node. The eigenvalues of the Jacobian matrix at (0, 0) are -1, 0 and this point is semi-hyperbolic. By using Theorem 2.19 of [5] we obtain that this point is a saddle. System (1) in the local chart (U2 , φ2 ) is u̇ = u2 + av 2 + u3 + (a − 1)uv 2 − buv 3 , v̇ = v(u2 + av 2 − bv 3 ). We see that the origin is an identically zero singular point of the above system. Applying the blow up (u, v) → (u, w), where u = u, v w = and rescaling the independent variable as ds = udt we obtain u the following system (3) u̇ = u(1+u+aw 2 +(a−1)uw 2 −bu2 w 3 ), ẇ = −w(1+aw 2 −uw 2 ). If a > 0, the origin is the unique singular point of system (3) on u = 0. The eigenvalues of the Jacobian matrix at this point are 1, -1 and the origin is a saddle. Going through the blow down we get the local behaviour at the origin of the chart U2 shown in Figure 4. If a < 0 system (1)√ has three singular points on u = 0, which are (0, 0) and p± = (0, ±1/ −a). The eigenvalues of the Jacobian matrix at (0, 0) are 1, -1 and it is a hyperbolic saddle. Computing the eigenvalues of the Jacobian matrix at the other singular points are 0, 2. These points are semi–hyperbolic and by using Theorem 2.19 of [5], we obtain that p± are unstable hyperbolic nodes if a < −b2 , p+ is a saddle-node and p− is an unstable node if a = −b2 , and p+ is a saddle and p− is an unstable node if a ∈ (−b2 , 0). Doing the blow down we obtain the local behaviour shown in Figures 5–7. 10 J. LLIBRE AND A. NABAVI v u 0 v w w a> u u u 0 > b a= 0 v b= 0 u Figure 4. Blowing down of the origin of the local chart U2 for a > 0. w w v u u u Figure 5. Blowing down of the origin of the local chart U2 for a < −b2 , b > 0. From the above information on the local phase portraits of the origin of U2 , we obtain the curves and the regions shown in Figure 8 with different behaviours. Actually, the origin of U2 is formed by • two hyperbolic sectors (topological index 0) in the region I, • two hyperbolic and two parabolic sectors (topological index 0) on the line {a = 0, b > 0}, SELKOV MODEL 11 w w v u u u Figure 6. Blowing down of the origin of the local chart U2 for a = −b2 , b 6= 0. w w v u u u Figure 7. Blowing down of the origin of the local chart U2 for a ∈ (−b2 , 0). • two elliptic and two parabolic sectors (topological index 2) in the region II, • one elliptic, two parabolic and one hyperbolic sectors (topolog- ical index 1) on the curve γ, • two hyperbolic and two parabolic sectors (topological index 0) in the region III, 3.3. On the periodic orbits. The curve β separates the region (3) from the region (4), see these regions in Figure 3. In the region (3) the 12 J. LLIBRE AND A. NABAVI a I b -0.25 γ III II Figure 8. The regions I, II, III, the point (b, a) = (0, 0), and the curve γ : {a = −b2 , b 6= 0} separating the different local phase portraits at the origin of the local chart U2 in the half plane {(b, a) : b ≥ 0, a ∈ R}. finite singular point is a hyperbolic stable focus, while in the region (4) it is a hyperbolic unstable focus. On the curve β and for b ∈ (0, b2 ] the finite singular point is a stable weak focus, and for b > b2 it is an unstable weak focus. So a limit cycle is born or disappears when the parameters (b, a) cross this curve β. In fact, crossing the curve β from the region (3) to (4) when b ∈ (0, b2 ) a stable limit cycle is born from the finite stable weak focus, while crossing the curve β from the region (4) to (3) when b > b2 an unstable limit cycle is born from the finite unstable weak focus. So the curve β is a curve of a Hopf bifurcation, for additional information on the Hopf bifurcations see for instance Chapter 3 of [6]. Numerically we have observed that these limit cycles which born in the described Hopf bifurcation end in a graphic having a part of it at infinity. Consequently, we have done the conjecture which appears in the section 1. The curves δ1 and δ2 in the boundaries of the regions R1 and R2 of Figure 2 have been obtained numerically, and for the values of the parameters on them the differential system (1) has a graphic where disappeared the limit cycle. Some approximations of the points (b, a) which lie on the curves δ1 and δ2 are given in Table 2. For any point (b, a) in the region R1 (shown in Figure 2) with 0 < b < b1 and a very close to zero, we can find numerically a stable periodic orbit. This periodic orbit does not exist when {a = 0, 0 < b < b1 }. SELKOV MODEL 13 a b on the curve 0 (0.9, 0.91) := b1 δ1 −0.001 (0.904, 0.905) δ1 −0.15 (1.168, 1.169) δ1 −1 (1.749, 1.75) δ2 −2 (2.154, 2.155) δ2 −3 (2.467, 2.468) δ2 Table 2. Some approximation of the points on the curves δ1 and δ2 a R3 P1 L3 R12 L1L R11 P3 b 2 L6 L4 P2 P4 L5 L8 R4 R5 L9 R6 R2 L7 Figure 9. Bifurcation diagram. On the other hand, since all the finite and infinite singular points are 1 √ unstable when 0 < a < (−1−2b2 + 1 + 8b2 ), the Poincaré–Bendixon 2 Theorem (see for instance Corollary 1.30 of [5]) shows that the system 1 has a stable limit cycle in the Poincaré disc when 0 < a < (−1 − √ 2 2b2 + 1 + 8b2 ), and we can conclude that the part of the curve δ1 with 0 < b < b1 is not only close to the b axis, actually it is the line {a = 0, 0 < b < b1 }. From the numerical results on the periodic orbits together with the information of Figure 8, we observe that the different behaviours in the regions I and III of Figure 8 cause that the region R1 splits into the regions R11 , R12 and the line L3 (see Figure 9). Then we can restate the conjecture of section 1 with more details as follows: 14 J. LLIBRE AND A. NABAVI P1 R3 -L1 L2 -P2 -L5 -L7 -R5 R11 L3 − R12 P3 -L6 -P4 -L8 -R4 -L4 R2 L9 R6 Figure 10. Local phase portraits in the Poincaré disc Conjecture. The differential system (1) has a unique stable limit cycle in the regions R11 , R12 and on the line L3 of Figure 9, and a unique unstable limit cycle in the region R2 of Figure 9. Moreover these limit cycles appear from a Hopf bifurcation at the point (b, b/(a + b2 )) when (b, a) are on the curve β and end in a a graphic with a part at infinity when (b, a) are on the curves δ1 and δ2 . 3.4. The Local phase portraits of all finite and infinite singular points in the Poincaré disc. Taking into account all the informa- tion of subsections 3.1, 3.2 and 3.3, we obtain the bifurcation diagram (shown in Figure 9), and we can draw in the Poincaré disc all the lo- cal phase portraits of the finite and infinite singular points (shown in Figure 10). SELKOV MODEL 15 Figure 11. The phase portraits of system (4). 3.5. The phase portraits in the Poincaré disc. At the point P1 , where (b, a) = (0, 0) we have the line of singular points x = 0 (and so an infinite number of separatrices), we can not use Theorem 2 for obtaining the global phase portrait. By rescaling ds = xdt we change system (1) (with a = b = 0) into the system (4) ẋ = −1 + xy, ẏ = −xy. Then easily we show that this system has the phase portrait shown in Figure 11. Adding the line of singular points x = 0 and changing the direction of the orbits in the region x < 0, we obtain the phase portrait of system (1) with a = b = 0 shown in Figure 1.A. By investigating all the local phase portraits in the Poincaré disc shown in Figure 10, and taking into account the Bendixson–Poincaré Theorem which states what are the α– and ω–limits of an orbit con- tained in a compact (in this case the compact is the Poincaré disc), it follows easily that the phase portraits of the differential system (1) are topologically equivalent to the phase portraits of Figure 1: • B when (b, a) ∈ R6 , • C when (b, a) ∈ L9 , • D when (b, a) ∈ {L2 , P2 , R5 , L5 , L7 }, • E when (b, a) ∈ {L3 , R12 }, • F when, (b, a) ∈ {P4 , L6 , L4 , R4 , P3 , L8 }, • G when (b, a) ∈ R2 , • H when (b, a) ∈ {R3 , L1 }, • I when (b, a) ∈ R11 . Acknowledgements The first author is partially supported by the Ministerio de Economı̀a, Industria y competitividad, Agencia Estatal de Investigación grant MTM2016-77278-P (FEDER), the Agència de Gestió d’Ajusts Univer- sitaris i de Recerca grant 2017 SGR 1617, and the European project 16 J. LLIBRE AND A. NABAVI Dynamics-H2020-MSCA-RISE-2017-777911. The second authoris sup- ported by Isfahan University of Technology (IUT). References [1] MJ. Álvarez, A. Ferragut and X. Jarque, A survey on the blow up technique, Int. J. Bifurc Chaos. 21 (2011), 3103–3118. [2] J.C. Artés, J. Llibre and C. Valls, Dynamics of the Higgings–Selkov and Selkov systems, Chaos, Solitons and Fractals. 114 (2018), 145–150. [3] P. Brechmann and A.D. Rendall, Dynamics of the Selkov oscillator, Math. Biosci. 306 (2018), 152–159. [4] H. Chen and Y. Tang, Proof of Artés–Llibre-Valls’s conjectures for the Higgins-Selkov systems, J. Differential Equations 266 (2019), 7638–7657. [5] F. Dumortier, J. Llibre and J.C. Artés, Qualitative Theory of Planar Differential Systems, Springer Verlag, New York, 2006. [6] Y. Kuznetsov, Elements of applied bifurcation theory, Applied Mathematical Sciences, Vol. 112, Springer-Verlag, New York, 2004. [7] L. Markus, Global structure of ordinary differential equations in the plane, Trans. Amer. Math. Soc . 76 (1954), 127–148. [8] D.A. Neumann, Classification of continuous flows on 2-manifolds, Proc. Amer. Math. Soc. 48 (1975), 73–81. [9] M.M. Peixoto, Proccedings of a simposium held at the university of Bahia, 349–420, Acad. Press, New York, 1973. [10] E. Selkov, Model of glycolytic oscillations, Eur. J. Biochem. 4 (1968), 79–86. [11] S.H. Strogatz, Nonlinear dynamics and chaos, Addison Wesley, 1994, pp. 205–207. 1 Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain E-mail address: jllibre@mat.uab.cat 2 Department of Mathematical Sciences, Isfahan University of Tech- nology, Isfahan, Iran, 84156-83111. E-mail address: a.nabavi@math.iut.ac.ir