EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6752 ISSN 1307-5543 – ejpam.com Published by New York Business Global Center Bifurcation for the Smallest Bimolecular Mass-Action System Rizgar H. Salih Department of Mathematics, College of Basic Education, University of Raparin, Rania, Kurdistan Region, Iraq Abstract. This paper investigates the center bifurcation of the smallest bimolecular mass-action system. A three-dimensional reaction network, consisting of three species and four reactions, governed by mass-action kinetics with a positive equilibrium point, is considered. In addition to the stability analysis of the equilibrium point, the dynamic directions of the model in its planes are examined. It has been previously shown that the equilibrium point is classified as a center when the reaction rate constants satisfy a specific condition, leading to a vertical Andronov-Hopf bifurcation. Furthermore, it is demonstrated that only one limit cycle can bifurcate from the center equilibrium point using a bifurcation technique. 2020 Mathematics Subject Classifications: 37G10, 37G15, 34C23 Key Words and Phrases: Bimolecular mass-action system, Hopf bifurcation, limit cycle, center bifurcation 1. Introduction Chemical reaction networks are fascinating for both practical and theoretical reasons. They are essential components of various biological models and significantly impact other fields of science and engineering. Numerous important findings, both traditional and con- temporary, provide insights into the dynamics of a chemical reaction network based on its combinatorial structure [1]. Wilhelm described a bimolecular chemical reaction network with three species and four reactions, having rank three, that exhibits Hopf bifurcation due to mass action dynamics [2]. His example demonstrated that the set of bimolecular networks with Hopf bifurcation is not empty. However, there are, up to isomorphism, 14670 bimolecular networks (with three species and four reactions of rank three) that permit positive equilibria [3]. The main question: is how many of these networks allow for Hopf bifurcation due to mass action dynamics? To answer this question, a scientific study by Banaji and Boros investigated the smallest DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6752 Email address: rizgar.salih@uor.edu.krd (R. H. Salih) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 2 of 14 bimolecular chemical reaction networks with Hopf bifurcation [3]. This means that these networks have the fewest species and reactions. The results show that the smallest bi- molecular mass-action network allowing Hopf bifurcation must have at least three species and at least four reactions. They found that 138 non-isomorphic bimolecular networks allow for Andronov-Hopf bifurcation. These networks can be dynamically divided into 87 non-equivalent classes. Of these classes, 86 admit non-degenerate Andronov-Hopf bifur- cations, leading to isolated limit cycles. However, in the remaining class, Andronov-Hopf bifurcations can only be degenerate. The three-dimensional system corresponding to the last class of the reaction network is ẋ1 = x1(k1x3 − k2x2), ẋ2 = x2(k2x1 − k3x3). (1) ẋ3 = −x3(k1x1 + k3x2) + 2k4, where ki, i = 1, 2, 3, 4 are real positive parameters called the reaction rate constants and xi ≥ 0. System (1) has a unique positive equilibrium point, which is of type center when the parameters satisfy k1 = k2+k3 [4]. It is proven by finding a constant of motion. They also proved the existence of a global center manifold that attracts all positive solutions. To gain information about bimolecular networks with Hopf bifurcation, the reader should consult the references [5–10]. An important question arises: If we perturb the parameters in system (1), how many periodic orbits can bifurcate from the positive equilibrium point? To address this, a technique is applied to system (1) to estimate the cyclicity bifurcating from the center. This method was first employed in three dimensions by Salih [11]. It has since been utilized in various scientific studies. Salih demonstrated that four limit cycles can bifurcate from the center of the 3D Lotka-Volterra system [11]. Additionally, Salih et al. [12] applied the technique to two differential systems, showing that one and five limit cycles can bifurcate from the quadratic polynomial system and the Lü system, respectively. Salih et al. [13] examined a specific type of the Jerk system, revealing that three and four limit cycles can bifurcate from the center equilibrium point under two different sets of conditions. This paper is organized as follows: Section 2 focuses on identifying the existence of equi- librium points in system (1) and studying their stability. The trajectory directions in the planes are investigated in Section 3. The next section is dedicated to the study of both Hopf and center bifurcations. Lastly, the conclusions are presented. 2. Equilibrium Point and Its Stability Equilibrium points are fundamental to understanding various aspects of dynamical systems. To identify the equilibrium points of system (1), the right-hand sides of the equations are set to zero. It is found that there is one positive equilibrium point, given by E = (√ k3k4 k1k2 , √ k1k4 k2k3 , √ k2k4 k1k3 ) . Through simple analysis, one can easily derive the following conclusion. R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 3 of 14 Proposition 1. For system (1), the stability of the equilibrium point E is determined as follows: i. It is asymptotically stable if and only if k1 > k2 + k3. ii. It is a center when k1 = k2 + k3. iii. It is unstable when k1 < k2 + k3. Proof. The Jacobian matrix of system (1) at the equilibrium points E is given by  0 − √ k2k3k4 k1 √ k1k3k4 k2√ k1k2k4 k3 0 − √ k1k3k4 k2 − √ k1k2k4 k3 − √ k2k3k4 k1 −2 √ k1k3k4 k2  (2) its characteristic equation is λ3 − Tλ2 −Kλ−D = 0, (3) where • T is the trace of (2) and T = −2 √ k1k3k4 k2 , • K is the sum of the diagonal minors of (2) and K = k4 (k3 − k1 − k2), • D is the determinant of (2) and D = −4k4 √ k1k2k3k4. Since D is negative, Eq. (3) has one negative real root and two other roots that have the same sign. As T is negative, the sign of TK+D is determined by (k1−k2−k3). The sign of (k1 − k2 − k3) plays an important role in determining the stability of the equilibrium point. Thus, i. When k1 > k2 + k3, it implies that TK +D > 0. According to the Routh-Hurwitz criterion, the real parts of the eigenvalues are negative. Therefore, the equilibrium point E is asymptotically stable if and only if k1 > k2 + k3, see Fig.1-a. ii. When k1 = k2+k3, Eq. (3) has one real negative root and a pair of purely imaginary roots. In reference [4], it has been proven that the equilibrium point E is of center type, see Fig.1-b. iii. When k1 < k2 + k3, it implies that TK +D < 0. In this case, Eq. (3) has two roots with positive real parts, making the equilibrium point E unstable, see Fig.1-c. R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 4 of 14 (a) (b) (c) Figure 1: The phase portrait of system (1) is shown and the parameters k2 = 1 2 , k3 = k4 = 1 are fixed. (a) when k1 = 2, E is asymptotically stable. (b) when k1 = 3 2 , E is a center. (c) when k1 = 1, E is unstable. The green and red balls indicate the equilibrium and initial points, respectively. 3. Trajectories in Planes This section focuses on the study of the dynamics of system (1) in the planes. To un- derstand the dynamics of the system, the trajectories and their evolution on the planes are investigated individually. Two of the planes, namely the x1x3−plane and the x2x3−plane, are invariant, meaning that the trajectories remain confined to them as time approaches infinity. Moreover, along all three coordinate axes, ẋ3 is always positive. On the x2x3−plane, ẋ2 is always negative. The horizontal isocline, x2x3 = 2k4 k3 , is used to determine the sign of ẋ3. A point (x2, x3) is chosen on one side of the isocline. First, a point on the left side is selected and x2x3 = 2k4 k3 − ϵ, where ϵ > 0 is set. This implies that ẋ3 = k3ϵ > 0. However, if a point on the other side is chosen, x2x3 = 2k4 k3 + ϵ, ϵ > 0, then ẋ3 = −k3ϵ < 0. The x3−axis is always invariant and on the x2−axis, ẋ3 > 0. Now, sufficient information has been gathered to sketch the trajectory directions of the system, as depicted in Fig. 2-a. On the x1x3-plane, ẋ1 is always positive and a horizontal isocline, x1x3 = 2k4 k1 , is present, which determines the sign of ẋ3. It can be easily shown that ẋ3 > 0 on the left side, while ẋ3 < 0 on the other side of the isocline. The x3-axis is found to be invariant and ẋ3 is positive on the x1-axis. By combining these observations, sufficient conditions for sketching the trajectory directions are derived and illustrated in Fig. 2-b. On the x1x2-plane, it is noted that ẋ1 is always negative, while both ẋ2 and ẋ3 are positive. The relations ẋ1+ ẋ2 = 0 and ẋ3 = 2k4 are satisfied, which implies that x1+x2 = c, where c ∈ R, and x3(t) = 2k4t. Since the plane is not invariant, this information indicates that the trajectories move with a positive slope. However, on each x1 and x2-axis, ẋ1 = ẋ2 = 0 and ẋ3 > 0. This indicates that near the axes, the trajectories move upward toward the x3-axis. This information helps us predict the direction of the trajectories. R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 5 of 14 By integrating the information obtained from Figs. 2-a, 2-b and the x1x2-plane, along with the dynamics at infinity, the complete information in Fig. 3 is derived. This infor- mation is useful for understanding the trajectory directions around the positive critical point E, which can be stable, unstable or a center in different cases. (a) (b) Figure 2: The trajectory directions of system (1) on the planes are shown. (a) trajectory direction on x2x3−plane. (b) trajectory direction on x1x3−plane. The dotted lines indicate the isoclines. Figure 3: Directions of the trajectory of system (1) in the three planes. 4. Bifurcation Analysis Bifurcations describe significant changes in the behavior of solution curves within a dynamical system as certain parameter values, called bifurcation values, are varied. This section examines both Hopf and center bifurcations in system (1), emphasizing the condi- tions under which they occur. R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 6 of 14 4.1. Hopf Bifurcation Consider the following dynamical system in three dimensions: Ẋ = f(X,µ), (4) where X is an element of R3, f is a real analytic function and µ belongs to Rk, serving as the bifurcation parameter. The essential conditions for the occurrence of a Hopf bifurcation are now presented. It is assumed that the system possesses an equilibrium point (E0, µ0) at which the following criteria are satisfied: i. The Jacobian matrix J = Df(E0, µ0) possesses a unique pair of purely imaginary eigenvalues λ(µ0) and λ̄(µ0), while the remaining eigenvalues are non-zero. ii. dRe(λ(µ0)) dµ ̸= 0. Subsequently, the system defined by (4) experiences a Hopf bifurcation at the equilibrium point (E0, µ0) [14]. System (4), satisfying the conditions above, can be transformed into the following canonical form: ẋ1 = −ωx2 + f1(x1, x2, x3;µ), ẋ2 = ωx1 + f2(x1, x2, x3;µ), (5) ẋ3 = λx3 + f3(x1, x2, x3;µ), where ω > 0, λ ̸= 0 and fi(x1, x2, x3;µ) = ∑∞ k=2 f k i (x1, x2, x3;µ) for i = 1, 2, 3 and fk i (x1, x2, x3;µ) are homogeneous polynomials of degree k. Proposition 2. For system (1), the Hopf bifurcation occurs at the equilibrium point E when the parameter k1 passes through k∗1 where k∗1 = k2 + k3. Proof. By applying the linear transformation x1 → x1 + √ k3k4 k1k2 , x2 → x2 + √ k1k4 k2k3 and x3 → x3 + √ k2k4 k1k3 , the equilibrium point E is relocated to the origin, resulting in the transformation of system (1) into: ẋ1 = (√ k3k4 k1k2 + x1 ) (−k2x2 + k1x3) , ẋ2 = (√ k1k4 k2k3 + x2 ) (k2x1 − k3x3) , (6) ẋ3 = − √ k1k2k4 k3 x1 − √ k2k3k4 k1 x2 − 2 √ k1k3k4 k2 x3 − k1x1x3 − k3x2x3. The characteristic equation corresponding to the Jacobian matrix of system (6) evaluated at the origin is λ3 + 2 √ k1k3k4 k2 λ2 − k4 (k3 − k1 − k2)λ+ 4k4 √ k1k2k3k4 = 0, (7) R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 7 of 14 Letting k1 = k∗1, then Eq. (7) can be rewritten as( λ+ 2 ω √ k3k4 (ω2 + 2k4k3) )( λ2 + ω2 ) = 0, where ω = √ 2k2k4. Clearly, Eq. (7) has two complex conjugate roots λ1,2 = ±iω and one real root λ3 = − 2 ω √ k3k4 (ω2 + 2k4k3). Thus, the initial requirement for the Hopf bifurcation theorem, which is the first condition, is fulfilled. It is important to note that, in general, λ = λ(k1). From Eq. (7), we can then express the relation as follows: f(λ(k1), k1) = ( λ(k1) + 2 ω √ k3k4 (ω2 + 2k4k3) )( λ2(k1) + ω2 ) . Then, the root λ = λ(k1) of Eq. (7) satisfies the following f(λ(k1), k1) = 0. (8) Differentiating Eq. (8) with respect to k1 yields ∂f ∂λ(k1) dλ(k1) dk1 + ∂f ∂k1 = 0, this implies that dλ(k1) dk1 = − ∂f ∂k1 ( ∂f ∂λ(k1) )−1 = − k2k4 ( λ2 + √ k1k2k4 k3 λ+ 2k2k4 ) √ k1k2k4 k3 ( 3k2λ2 + 4 √ k1k2k3k4 λ+ k2k4 (k1 + k2 − k3) ) . (9) By taking the root λ(k∗1) = λ1,2(k ∗ 1), we find that λ1,2(k ∗ 1) = ±iω. Substituting the value of λ1 into Eq. (9), we obtain: d dk1 Re(λ1,2(k1))|k1=k∗1 = − k4ω √ k3k4 (ω2 + 2k3k4) ω4 + 4k3k4 (ω2 + 2k3k4) ̸= 0. Thus, the second condition for a Hopf bifurcation is met. Consequently, at E, system (1) undergoes a Hopf bifurcation when k1 = k∗1. It is straightforward to determine the Lyapunov coefficients corresponding to system (6) at k1 = k∗1. Since the equilibrium point is a center, the k−th Lyapunov coefficients vanish for all k at k∗1 [15]. Therefore, by the Lyapunov theorem, system (1) exhibits a vertical Andronov-Hopf bifurcation as k1 passes through k∗1. 4.2. Center Bifurcation Research in bifurcation theory currently focuses on the bifurcation of limit cycles from critical points. A limit cycle can be achieved by perturbing a focus or center. A widely used method is center bifurcation, which helps estimate cyclicity and examine the bifurcation of limit cycles from the center (see [16] and [17]). Here, we consider system (5). The set of all parameters in fi(x1, x2, x3) for i = 1, 2, 3 is denoted by Λ and K is the corresponding parameter space. R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 8 of 14 4.2.1. A Technique to Examine the Cyclicity This subsection describes a technique used for analyzing cyclicity. As previously noted, Christopher [18] explored a method for analyzing cyclicity bifurcating from a center in two-dimensional systems by linearizing the Lyapunov quantities. Salih [11] extended this approach to three-dimensional systems to estimate the cyclicity of the center, applying it for the first time to three-dimensional Lotka-Volterra systems. The technique used to estimate cyclicity in three-dimensional differential systems can be summarized in the following steps: (i) Select a point on a center variety. (ii) Linearize the Lyapunov quantities around this point. (iii) Determine the codimension of the point. If the codimension of the selected point on the center variety is r and the first r linear terms of the Lyapunov quantities are linearly independent, then the cyclicity is r − 1. This means that r − 1 limit cycles can bifurcate from a small perturbation. Defining the Lyapunov function and computing its focal values is a traditional approach to assess the number of limit cycles and their stability. This method involves looking for a function of the following form: F (x1, x2, x3) = x21 + x22 + ∞∑ k=3 Fk(x1, x2, x3), (10) where Fk = ∑k i=0 ∑i j=0Ck−i,i−j,jx k−i 1 xi−j 2 xj3 for system (5) and the coefficients of Fk sat- isfy X (F ) = L1(x 2 1 + x22) + L2(x 2 1 + x22) 2 + L3(x 2 1 + x22) 3 + ..., (11) where Li, i = 1, 2, .. are polynomials in the parameters of the system and the Li is called the i−th Lyapunov constant (focal value). To explain the technique in greater detail, it is assumed that the center critical point of (5) corresponds to 0 ∈ K, using a perturbation method in the parameters. This can be written: X = X0 + X1 + X2 + ..., F = F0 + F1 + F2 + ..., (12) Li = Li0 + Li1 + Li2 + ..., i = 1, 2, ..., where X0, F0 and Li0 are calculated at the unperturbed parameters, X1, F1 and L1i are obtained at a perturbed parameters of first order (they contain the terms of degree one in Λ), X2, F2 and L2i are obtained at a perturbed parameters of second order (they contain the terms of degree two in Λ) and so forth. The Lyapunov function Fi and the Lyapunov R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 9 of 14 quantity Li are both of degree i in terms of parameters. By substituting Eq. (12) into Eq.(11), we obtain: X0F0 = 0, X0F1 + X1F0 = L11(x 2 1 + x22) + L21(x 2 1 + x22) 2 + ... , X0F2 + X1F1 + X2F0 = L12(x 2 1 + x22) + L22(x 2 1 + x22) 2 + ... (13) and more general, X0Fi + ...+ XiF0 = L1i(x 2 1 + x22) + L2i(x 2 1 + x22) 2 + ... , i = 1, 2, 3, ... (14) The linear terms of the Lyapunov quantities Lk (modulo the Li, i < k) can be obtained by simultaneously solving the two equations in (13) using linear algebra. Eq. (14) is then used to derive the higher-order terms of the Lyapunov quantities. 4.2.2. Cyclicities Bifurcated from the Center Equilibrium Point To apply the above technique, we choose point (k1, k2, k3, k4) = (32 , 1 2 , 1, 1) on center variety and let k1 = 3 2 + a1 + a2, k2 = 1 2 + b1 + b2, k3 = 1 + c1 + c2 and (15) k4 = 1 + d1 + d2 are defined where a1, b1, c1, d1 and a2, b2, c2, d2 are parameters introduced by perturbation in the system of first and second order, respectively. Therefore, the unperturbed vector field X0, the first-order perturbed vector field X1 and the second-order perturbed vector field X2 are defined as follows: X0 = 1 6 ( 2 √ 3 + 3x1 ) (3x3 − x2) ∂ ∂x1 + 1 2 (x1 − 2x3) ( x2 + √ 3 ) ∂ ∂x2 − 1 6 ( 3 √ 3 x1 + 2 √ 3x2 + 12 √ 3x3 + 9x1x3 + 6x2x3 ) ∂ ∂x3 , X1 = (e2x2 + e3x3 − b1x1x2 + a1x1x3) ∂ ∂x1 + (e1x1 − e3x3 + b1x1x2 − c1x2x3) ∂ ∂x2 + (e2x2 − 2e3x3 − a1x1x3 − c1x2x3) ∂ ∂x3 , X2 = (f2x2 + f3x3 + a2x1x3 − b2x1x2) ∂ ∂x1 + (f1x1 − f3x3 + b2x1x2 − c2x2x3) ∂ ∂x2 + (f2x2 − 2f3x3 − a2x1x3 − c2x2x3) ∂ ∂x3 , (16) where e1 = √ 3 12 (2a1 + 6b1 − 3c1 + 3d1) , e2 = √ 3 18 (2a1 − 6b1 − 3c1 − 3d1) , e3 = √ 3 6 (2a1 − 6b1 + 3c1 + 3d1) , R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 10 of 14 f1 = √ 3 144 (−4a21 + 24b1a1 − 12a1c1 + 12a1d1 − 36b21 − 36b1c1 + 36b1d1 + 27c21 − 18c1d1 − 9d21 + 24a2 + 72b2 − 36c2 + 36d2), f2 = √ 3 72 (−4a21 + 8b1a1 + 4a1c1 + 4a1d1 + 12b21 − 12b1c1 − 12b1d1 + 3c21 − 6c1d1 + 3d21 + 8a2 − 24b2 − 12c2 − 12d2) and f3 = √ 3 72 (−4a21 − 24b1a1 + 12a1c1 + 12a1d1 + 108b21c− 36b1c1 − 36b1d1 − 9c21 + 18c1d1 − 9d21 + 24a2 − 72b2 + 36c2 + 36d2). The main result regarding the bifurcated limit cycle from the center equilibrium point, using the technique described above, is the following theorem. Theorem 1. For system (1), when k1 = k2 + k3, using first-order perturbation, only a single limit cycle can bifurcate from the equilibrium point at E. Proof. Instead of analysing system (1) at E, system (6) is examined at the origin. Using the linear transformation X = PY, P = −2 √ 3 −2 −1/6 2 √ 3 −3 1/4 0 1 5/12  , (17) where X = (x1, x2, x3), Y = (x, y, z), the linear part of system (6) at the origin can be written in the real canonical form as0 −1 0 1 0 0 0 0 −2 √ 3  , and the new system is given by ẋ = −y − √ 3x2 + xy + 7 √ 3 13 y2 + 5 √ 3 156 yz + √ 3 312 z2, ẏ = x+ √ 3xy + 5 √ 3 12 xz + 3 13 y2 + 5 26 yz + 25 624 z2, (18) ż = −2 √ 3z + 180 13 y2 + 72 13 yz − 5 52 z2. The transformation described in Eq. (17) is applied to the first order perturbed vector field component of system (6), leading to the following results: ẋ =− 1 13 (4e1 + 15e2)x− √ 3 78 (8e1 − 45e2 + e3)y − √ 3 936 (8e1 + 45e2 + 5e3)z + 1 13 (3a1 + 13b1 + 2c1)xy + 1 156 (15a1 − 65b1 + 10c1)xz + √ 3 78 (3a1 − 2c1)yz − 2 √ 3b1x 2 + √ 3 13 (a1 + 13b1 − c1)y 2 R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 11 of 14 + √ 3 1872 (5a1 − 13b1 + 5c1)z 2, ẏ = √ 3 13 (5e1 − 4e2)x+ 1 13 (5e1 + 6e2 − e3)y + 1 156 (5e1 − 6e2 − 5e3)z + 2 √ 3 13 (3a1 + 2c1)xy + 5 √ 3 78 (3a1 + 2c1)xz + 1 13 (3a1 − 2c1)yz + 6 13 (a1 − c1)y 2 + 5 312 (a1 + c1)z 2, ż = −12 √ 3 13 (e1 − 6e2)x− 1 13 (12e1 + 108e2 + 60e3)y − 1 13 (e1 − 9e2 + 25e3)z + 24 √ 3 13 (2a1 − 3c1)xy + 10 √ 3 13 (2a1 − 3c1)xz + 12 13 (2a1 + 3c1)yz + 12 13 (4a1 + 9c1)y 2 + 5 156 (4a1 − 9c1)z 2. (19) Now, the unperturbed Lyapunov function, F0, and the first-order perturbed Lyapunov function, F1, are defined by: F0 = x2 + y2 + N∑ k=3 k∑ i=0 i∑ j=0 Ck−i,i−j,jx k−iyi−jzj , F1 = N∑ k=3 k∑ i=0 i∑ j=0 Dk−i,i−j,jx k−iyi−jzj , (20) where N ≥ 3. It is easy to show that the Lyapunov function, F0, of Eq. (18) is satis- fied by X0F0 = 0. Using the computer algebra package MAPLE, the following linearly independent terms of Lyapunov quantities are given by Eq. (13): (i) L11 = − √ 3 156(14a1 − 54b1 − 9c1 − 15d1). (ii) L21 = 3 √ 3 21632(2254a1 − 4534b1 − 1969c1 − 855d1). The origin of system (6) is weak focus of order one if and only if a1 = 1 14 (54b1 + 9c1 + 15d1). (21) Since the Jacobian of L11 and L21 with respect to a1 and b1 is non-zero, it indicates that, by suitable perturbation of the coefficients of the Lyapunov quantities, only one limit cycle can be bifurcated from the equilibrium point E of system (1) in the neighborhood of that point. If we take b1 = c1 = d1 = 0.1 as a numerical example, then from equation (21), we obtain a1 = 39 70 and a limit cycle is observed from the first-order perturbation. From Fig. 4, we note that the blue trajectory moves inward toward the equilibrium point, while the red trajectory moves outward. This predicts that there may be a limit cycle located in a region between the two initial points, although determining the existence of the limit cycle is not an easy task. R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 12 of 14 Figure 4: The phase portrait of system (1) around the positive equilibrium point is shown. The red and green balls indicate the equilibrium point and the initial points, respectively. It is also observed that applying the second-order perturbation in Eq. (14) does not change the outcome regarding the number of perturbed limit cycles, which remains one. Understanding the uniqueness of the limit cycle in bimolecular mass-action systems has important implications for chemical reaction dynamics and computational chemistry. This framework enhances our comprehension of chemical behavior under varying conditions, aiding in the prediction of real-world reaction outcomes. Additionally, these findings can improve computational models, enhancing the accuracy of simulations for complex chemical networks. By connecting theory with practical applications, this research fosters a deeper understanding of chemical systems and supports future advancements in the field. 5. Conclusions In this paper, the smallest bimolecular mass-action system exhibiting a center equi- librium point is studied. It has previously been shown that when a key reaction rate parameter k1 satisfies k1 = k2 + k3, the positive equilibrium point was classified as a center and a vertical Andronov-Hopf bifurcation was observed. The novelty of this re- search lies in determining how many limit cycles can bifurcate from the center equilibrium point. Using a perturbation technique, the number of limit cycles bifurcating from the center equilibrium point is estimated and it is concluded that only one limit cycle can emerge. These findings enhance the understanding of the dynamics of bimolecular reac- tion networks and provide a foundation for further research on the bifurcation behavior of R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 13 of 14 complex chemical systems. The presence of the limit cycle is important because it offers valuable insights into the stability and behavior of chemical reaction networks, forming a foundation for comprehending more complex systems. Future studies could focus on higher-dimensional mass-action systems, which might display more complex bifurcation patterns. Furthermore, examining the implications of these results in practical chemical processes and integrating them into computational models could deepen our understanding and application of dynamical systems in chemistry. References [1] Rachel S Lawrence. Simplicial Reaction Networks and Dynamics on Graphs. PhD thesis, University of California, Berkeley, 2023. [2] Thomas Wilhelm. The smallest chemical reaction system with bistability. BMC systems biology, 3:1–9, 2009. [3] Murad Banaji and Balázs Boros. The smallest bimolecular mass action reaction networks admitting andronov–hopf bifurcation. Nonlinearity, 36(2):1398, 2023. [4] Murad Banaji, Balázs Boros, and Josef Hofbauer. The smallest bimolecular mass- action system with a vertical andronov–hopf bifurcation. Applied Mathematics Let- ters, 143:108671, 2023. [5] Balázs Boros and Josef Hofbauer. Limit cycles in mass-conserving deficiency-one mass-action systems. Electron. J. Qual. Theory Differ. Equ., (42):1–18, 2022. [6] Péter Érdi and János Tóth. Mathematical models of chemical reactions: theory and applications of deterministic and stochastic models. Manchester University Press, 1989. [7] Gy Póta. Two-component bimolecular systems cannot have limit cycles: A complete proof. The Journal of Chemical Physics, 78(3):1621–1622, 1983. [8] György Póta. Irregular behaviour of kinetic equations in closed chemical systems. os- cillatory effects. Journal of the Chemical Society, Faraday Transactions 2: Molecular and Chemical Physics, 81(1):115–121, 1985. [9] Thomas Wilhelm and Reinhart Heinrich. Smallest chemical reaction system with hopf bifurcation. Journal of mathematical chemistry, 17(1):1–14, 1995. [10] Thomas Wilhelm and Reinhart Heinrich. Mathematical analysis of the smallest chemical reaction system with hopf bifurcation. Journal of Mathematical Chemistry, 19(2):111–130, 1996. [11] Rizgar Salih. Hopf bifurcation and centre bifurcation in three dimensional Lotka- Volterra systems. PhD thesis, Plymouth University, 2015. [12] Rizgar Salih and Mohammad Hasso. Centre bifurcations of periodic orbits for some special three dimensional systems. Electronic Journal of Qualitative Theory of Dif- ferential Equations, 2017(19):1–10, 2017. [13] Rizgar Salih, Mohammad Hasso, and Surma Ibrahim. Centre bifurcations for a three dimensional system with quadratic terms. Zanco Journal of Pure and Applied Sci- ences, 32(2):62–71, 2020. R. H. Salih / Eur. J. Pure Appl. Math, 18 (4) (2025), 6752 14 of 14 [14] John Guckenheimer and Philip Holmes. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, volume 42. Springer Science & Business Media, 2013. [15] Valery Romanovski and Douglas Shafer. The center and cyclicity problems: a com- putational algebra approach. Springer Science & Business Media, 2009. [16] Nikolai Nikolaevich Bautin. On the number of limit cycles appearing with varia- tion of the coefficients from an equilibrium state of the type of a focus or a center. Matematicheskii Sbornik, 72(1):181–196, 1952. [17] P Yu and M Han. Twelve limit cycles in a cubic order planar system with z˜ 2 symmetry. Communications on pure and applied analysis, 3:515–526, 2004. [18] Colin Christopher. Estimating limit cycle bifurcations from centers. In Differential equations with symbolic computation, pages 23–35. Springer, 2005.