Chemical Reaction Systems with a Homoclinic Bifurcation: an Inverse Problem

An inverse problem framework for constructing reaction systems with prescribed properties is presented. Kinetic transformations are defined and analysed as a part of the framework, allowing an arbitrary polynomial ordinary differential equation to be mapped to the one that can be represented as a reaction network. The framework is used for construction of specific two- and three-dimensional bistable reaction systems undergoing a supercritical homoclinic bifurcation, and the topology of their phase spaces is discussed.


Introduction
Chemical reaction networks under the mass action kinetics are relevant for both pure and applied mathematics. The time evolution of the concentrations of chemical species is described by kinetic equations which are a subset of first-order, autonomous, ordinary differential equations (ODEs) with polynomial right-hand sides (RHSs). On the one hand, the kinetic equations define a canonical form for a subset of ODEs, thus being important for pure mathematics [1,2]. They can display not only the chemically regular phenomenon of having a globally stable fixed point, but also the chemically exotic phenomena (multistability, oscillations and chaos). It is then no surprise that chemical reaction networks can perform the same computations as other types of physical networks, such as electronic and neural networks [3]. On the other hand, reaction networks are a versatile modelling tool, decomposing processes from applications into a set of simpler elementary steps (reactions). The exotic phenomena in systems biology often execute specific biological functions, example being the correspondence between limit cycles and biological clocks [4,5].
The construction of reaction networks displaying prescribed properties may be seen as an inverse problem in formal reaction kinetics [6], where, given a set of properties, a set of compatible reaction networks is searched for. Such constructions are useful in many application areas, including systems biology (as caricature models), synthetic biology (as blueprints), and numerical analysis (as test problems) [7][8][9]. In systems biology, kinetic ODEs often have higher nonlinearity degree and higher dimension, thus not being easily amenable to mathematical analysis. Having ODEs with lower nonlinearity degree and lower dimension allows for a more detailed mathematical analysis, and also adds to the set of test problems for numerical methods designed for more challenging real-world problems. In synthetic biology, such constructed systems may be used as a blueprint for engineering artificial networks [9,10].
A crucial property of the kinetic equations is a lack of so-called cross-negative terms [11], corresponding to processes that involve consumption of a species when its concentration is zero. Such terms are not directly describable by reactions, and may lead to negative values of concentrations. The above mentioned issues with cross-negative terms, together with the requirement that the dependent variables are uniformly bounded in time, imply that not every nonnegative polynomial ODE system is kinetic, and, thus, further constrain the possible dynamics. A trivial example of an ODE with a cross-negative term is given by dx/dt = −k, for constant k > 0, where the term −k, although a polynomial of degree zero, nevertheless cannot be directly represented by a reaction, and results in x < 0.
In two dimensions, where the phase plane diagram allows for an intuitive reasoning, the exotic dynamics of ODE systems reduces to cycles and multistability. While two-dimensional nonkinetic polynomial ODE systems exhibiting a variety of such dynamics can be easily found in the literature, the same is not true for the more con-strained two-dimensional kinetic ODE systems [12] (see also Appendix 1). Motivated by this, this paper consists of two main results: firstly, building upon the framework from [6,11], an inverse problem framework suitable for constructing the reaction systems is presented in Sect. 3, with the focus on the so-called kinetic transformations, allowing one to map a nonkinetic into a kinetic system. Secondly, in Sect. 4, the framework is used for construction of specific two-and three-dimensional bistable kinetic systems undergoing a global bifurcation known as a supercritical homoclinic bifurcation. The corresponding phase planes contain a stable limit cycle and a stable fixed point, with a parameter controlling the distance between them, and their topology is discussed. Definitions and basic results regarding reaction systems are presented in Sect. 2. A summary of the paper is presented in Sect. 5.

Notation and definitions
The notation and definitions in this paper are inspired by [11,13,14]. The formal sum notation is introduced so that unnecessary ordering of elements of a set can be avoided, such as when general frameworks involving sets are described, and when objects under consideration are vector components with irrelevant ordering. The usual vector notation is used when objects under consideration are equations in matrix form, and is put into using starting with Eq. (9).

Reaction networks and reaction systems
Definition 2.2 A reaction network is a triple {S, C, R}, where (i) S is a finite set, with elements s ∈ S called the species of the network.
(ii) C ⊂ N S is a finite set, with elements c ∈ C, called the complexes of the network, such that c∈C supp(c) = S. Components of c are called the stoichiometric coefficients. (iii) R ⊂ C × C is a binary relation with elements r = (c, c ), denoted r = c → c , with the following properties: (a) ∀c ∈ C (c → c) / ∈ R; Elements r = c → c are called reactions of the network, and it is said that c reacts to c , with c being called the reactant complex, and c the product complex.
The order of reaction r is given by o r = s∈S c s < ∞ for r = c → c ∈ R.
Note that as set R implies sets S and C, reaction networks are denoted with R, for brevity. Also, as it is unlikely that a reaction between more than three reactants occurs [11], in this paper we consider reactions with o r ≤ 3. To represent some of the non-chemical processes as quasireactions, the zero complex is introduced, denoted with ∅, with the property that supp(∅) = ∅, where ∅ is the empty set.

Definition 2.3
Let R be a reaction network and let κ : R S ≥ → R R ≥ be a continuous function which maps x ∈ R S ≥ (called "species concentrations") into κ(x) ∈ R R ≥ (called "reaction rates"). Then κ is said to be a kinetics for R provided that, for all x ∈ R S ≥ and for all r = (c → c ) ∈ R, positivity κ r (x) > 0 is satisfied if and only if supp(c) ⊂ supp(x).
An interpretation of Definition 2.3 is that a reaction, to which a kinetics can be associated, can occur if and only if all the reactant species concentrations are nonzero.

Definition 2.4
A reaction network R augmented with a kinetics κ is called a reaction system, and is denoted {R, κ}. Definition 2.5 Given a reaction system {R, κ}, the induced kinetic function, K(·; R) : The induced system of kinetic equations, describing the time evolution of species concentrations x ∈ R S ≥ , takes the form of a system of autonomous first-order ordinary differential equations (ODEs), and is given by Note that the kinetic function uniquely defines the system of kinetic equations, and vice-versa. In this paper, the species concentrations satisfying equation (1) are required to be finite, i.e. x s < ∞, for s ∈ S, and for t ≥ 0. Infinite-time blow-ups are allowed for initial conditions located on a finite number of (S − z)-dimensional subspaces of R S ≥ , z ≥ 1. Finite-time blow-ups [15], being chemically more unrealistic, are not allowed in this paper.

Definition 2.6
Kinetics κ is called the mass action kinetics if κ r (x) = k r x c , for r = (c → c ) ∈ R, where k r > 0 is the rate constant of reaction r , and x c = s∈S x c s s , with 0 0 = 1. A reaction system with the mass action kinetics is denoted {R, k}, and the corresponding kinetic function is denoted A review of the mass action kinetics can be found in [16]. In this paper, most of the results are stated with kinetics fixed to the mass action kinetics. Let us note that the results are also applicable to nonpolynomial ODEs that can be reduced to the polynomial ones [1].
Example 1 Consider the following reaction network (consisting of one reaction) under the mass action kinetics: so that S = {s 1 , s 2 }, C = {s 1 + s 2 , 2s 2 }, R = {s 1 + s 2 → 2s 2 } and k = {k 1 }. Concentration x ∈ R S ≥ has two components. To simplify our notation, we write x 1 = x s 1 , x 2 = x s 2 , and K 1 (x; k) = K s 1 (x; R), K 2 (x; k) = K s 2 (x; R). Then the induced system of kinetic equations is given by

Kinetic and nonkinetic functions
In this subsection, nonkinetic functions are defined, and further notation for kinetic and nonkinetic functions taking the mass action form is presented.
≥ , s ∈ S and r ∈ R. If ∃s ∈ S, ∃r ∈ R and ∃x ∈ R S ≥ such that s ∈ supp c (x) and f sr (x) < 0, then f sr (x) is called a cross-negative term, and function f (x) and ODE system dx/dt = f (x) are said to be nonkinetic.
An interpretation of a cross-negative term is that the process corresponding to such a term would consume at least one reactant even when its concentration is zero, so that it cannot be represented as kinetic reactions.
Kinetic and nonkinetic functions taking the mass action kinetics form are central to this paper. The related notation is introduced in the following definition. Definition 2.8 Let P(·; k) : R S → R S , k ∈ R R , be a polynomial function with polynomial degree deg(P(x; k)) ≤ m, m ∈ N, i.e. P s (x; k) is a polynomial in x ∈ R S of degree at most m with coefficients k ∈ R R , for all s ∈ S. Then, the set of functions P(x; k) is denoted by P m (R S ; R S ). If P(x; k) is a kinetic function, it is denoted by K(x; k), k ∈ R R > , and the set of such functions is denoted by P K m (R S ≥ ; R S ). If P(x; k) is a nonkinetic function, it is denoted by N (x; k), k ∈ R R , and the set of such functions with domain R S is denoted by P N m (R S ; R S ), while with domain R S ≥ by P N m (R S ≥ ; R S ). Note that a system {R, k}, corresponding to N (x; k) in Definition 2.8, has a welldefined reaction network R (for r = c → c , r ∈ R, we restrict c, c to positive integers), but an ill-defined kinetics taking the mass action form (we allow set k to have elements that are negative). Thus, set k corresponding to N (x; k) cannot be interpreted as a set of reaction rate constants, as opposed to set k corresponding to K(x; k) (see also Example 2). Note also that P m (R S ≥ ;

Properties of kinetic functions
From Definition 2.3 it follows that a kinetic function K(x; R) has a structural property: cross-negative terms are absent. In this subsection, further properties of K(x; R) are defined: nonnegativity (absence of cross-negative effect), and a structural property called x-factorability.
≥ 0, then f (x) and dx/dt = f (x) are said to be nonnegative. Otherwise, f (x) and dx/dt = f (x) are said to be negative, and a cross-negative effect is said to exists ∀x ∈ R S ≥ for which ∃s ∈ S such that s ∈ supp c (x) and f s (x) < 0.
Note that the absence of cross-negative terms implies nonnegativity, but the converse is not necessarily true [11,17], i.e. an ODE system may have cross-negative terms, without having a cross-negative effect, as we will show in Example 2. Cross-negative terms play an important role in mathematical constructions of reaction systems, in the context of chaos in kinetic equations, and pattern formation via Turing instabilities [18,19]. It has been proved that quadratic two-dimensional ODEs without cross-negative terms, and of the form such that the induced reactions involve at most two reactants and two products, cannot display limit cycle oscillations [12,20]. In Appendix 1, by generalizing [12], we prove that the nonexistence of a cross-negative effect in the ODEs of such form is a sufficient condition for nonexistence of limit cycles.
System (5)-(6) induces a reaction system only in case (i). In particular, nonnegative ODE system (5)-(6) with P(x; k) ∈ P N 2 (R S ≥ ; R S ) in part (ii)(a) does not induce a reaction system (although, given a nonnegative initial condition, the solution of (5)- (6) is nonnegative for all forward times).
≥ , s ∈ S and r ∈ R. Then term f sr (x) is said to be x s -factorable if f sr (x) = x s p sr (x), where p sr (x) is a polynomial function of x. If ∃s ∈ S, such that f s (x) = x s r ∈R p sr (x), then f (x) and ODE system dx/dt = f (x) are said to be x s -factorable. If ∀s ∈ S it is true that f s (x) = x s r ∈R p sr (x), then f (x) and ODE system dx/dt = f (x) are said to be x-factorable.
X -factorable ODE systems are a subset of kinetic equations under the mass action kinetics [21] (see also Sect. 3.2.2).

Inverse problem for reaction systems
In some applications, we are interested in the direct problem: we are given a reaction network with kinetics, i.e. a reaction system {R, κ}, and we then analyse the induced system of kinetic equations (1) in order to determine properties of the reaction system. For example, an output of a direct problem might consist of verifying that the kinetic equations undergo a bifurcation. In this paper, we are interested in the inverse problem: we are given a property of an unknown reaction system, and we would then like to construct a reaction system displaying the property. The inverse problem framework described in this section is inspired by [6,11].
The first step in the inverse problem is, given a quantity that depends on a kinetic function, to find a compatible kinetic function K(x; R), while the second step is then to find a reaction system {R, κ} induced by the kinetic function. The second step is discussed in more detail in Sect. 3.1, while the first step in Sect. 3.2. The constructions of a reaction system {R, κ} often involve constraints defining simplicity of the system (e.g. see [22]), and the simplicity can be related to the kinetic equations (structure and dimension of the equations, and/or the phase space), and/or to reaction networks (cardinality, conservability, reversibility, deficiency). How the simplicity constraints are prioritized depends on the application area, with simplicity of the kinetic equations being more important for mathematical analysis, while simplicity of the reaction networks for synthetic biology.

The canonical reaction network
Let us assume that we are able to construct an ODE system of the form (1) where its RHS is a kinetic function, K(x; R), and the system has the property required by the inverse problem. Then, one can always find a reaction system induced by the kinetic function [6,18]. While, for a fixed kinetics, a reaction network induces a kinetic function uniquely by definition (see Definition 2.5), the converse is not true-the inverse mapping of the kinetic function to the reaction networks is not unique-a fact known as the fundamental dogma of chemical kinetics [14,18,23,24]. For example, in [24], for a fixed kinetic function and a fixed set of complexes (C fixed), mixed integer programming is used for numerical computation of different induced reaction networks with varying properties. On the other hand, a constructive proof that every kinetic function induces a reaction system is given in [6,18], where C is generally not fixed (product complexes may be created), but the construction can be performed analytically, and it uniquely defines an induced reaction system for a given kinetic function under the mass action kinetics. The same procedure is used in this paper, so it is now defined.
Consider the kinetic function given by ≥ and d r ∈ R S . Let us map K(x; R) to a reaction system {R K −1 , κ K −1 } with complexes and kinetics given by: is not allowed to induce reactant complexes. Then, the complexes are uniquely obtainable in the case of the mass action kinetics. (ii) Reaction c r → c rs is then constructed for each r ∈ R and s ∈ S, where new product complexes are given by c rs = c r + sign(d rs )s, with sign(·) being the sign function.
Note that the procedure in Definition 3.1 creates a reaction for each term in each kinetic equation. Note also that each reaction leads to a change in copy number of precisely one chemical species, and the change in the copy number is equal to one. Thus, the canonical reaction networks are simple in the sense that they can be constructed from a kinetic function in a straightforward way, while they generally do not contain minimal number of reactions, nor complexes.

Example 4
The canonical reaction network for system (3)-(4) is given by Note that the canonical reaction network (7) contains more reactions than the original network (2).

Kinetic transformations
Firstly, mapping a solution-dependent quantity to the RHS of an ODE system is much more likely to result in nonkinetic functions, N (x; R), on the RHS (see Definition 2.7) [11]. However, only kinetic functions induce reaction networks, as exemplified in Example 2. Secondly, even if mapping a solution-dependent quantity results in a kinetic function, it may be necessary to modify the function in order to satisfy given constraints, and this may change the kinetic function into a nonkinetic function. For these two reasons, it is beneficial to study mappings that can transform arbitrary functions into kinetic functions. This motivates the following definition, for the case of mass action kinetics, that relies on the notation introduced in Definition 2.8.
Consider the corresponding ODE system in the formal sum notation where x ≡ x(t) ∈ R S . Then, a transformation Ψ is called a kinetic transformation if the following conditions are satisfied: (ii) Let x * be a fixed point of (8) that is mapped by Ψ to fixed pointx * ∈ RS ≥ of the system of kinetic equations (1) with K(x;k) on its RHS. Let also the eigenvalue of the Jacobian matrix of P(x; k), J (x * ; k), denoted by λ n , be mapped to the eigenvalue of Jacobian of K(x;k), J Ψ (x * ;k), which is denoted byλ n , for n = 1, 2, . . . , S. Then, for every such fixed point x * , it must be true that sign(Re(λ n )) = sign(Re(λ n )) (preservation of the fixed point stability), sign(Im(λ n )) = sign(Im(λ n )) (preservation of the fixed point orientability), n = 1, 2, . . . , S, and, if there are any additional eigenvaluesλ n , n = S + 1, S + 2, . . . ,S, they must satisfy sign(Re(λ n )) < 0 (stability of the fixed point along the added dimensions), where Re(·) and Im(·) denote the real and imaginary part of a function.
If any of the condition (i)-(ii) is not true, Ψ is called a nonkinetic transformation.
Put more simply, given an input polynomial function, a kinetic transformation must (i) map the input polynomial function into an output kinetic function, and (ii) the output function must be linearly locally topologically equivalent to the input function in the neighbourhood of the corresponding fixed points, and the fixed points are asymptotically stable along any additional dimensions of the output function (corresponding to the additional species), preventing the state of the system from escaping the fixed point neighbourhood along the added dimensions. Let us note that the output function is defined only in the nonnegative orthant, so that the topological equivalence must hold only near the fixed points of the input function that are mapped to the nonnegative orthant under kinetic transformations. One may wish to impose a set of constraints on an output function, such as requiring that a predefined region of interest in the phase space of the input function is mapped to the positive orthant of the corresponding output function. The subset of constraints is now defined that can be formulated purely in terms of the reaction rate constants.
There are two sets of kinetic transformations. The first, and the preferred, set of possible kinetic transformations are affine transformations, which are discussed in Sect. 3.2.1. Affine transformations may be used, not only as possible kinetic transformations, but also to satisfy a set of constraints. The second set, necessarily used when affine transformations fail, are nonlinear transformations that replace cross-negative terms, with x-factorable terms (see Definition 2.10), without introducing new crossnegative terms, and two such transformations are discussed in Sects. 3.2.2 and 3.2.3. In choosing a nonlinear transformation, one generally chooses between obtaining, on the one hand, lower-dimensional kinetic functions with higher-degree of nonlinearity [i.e. lowerS/S and higherm/m in Definition 3.2(i)] and/or higher numbers of the nonlinear terms, and, on the other hand, higher-dimensional kinetic functions with lower degree of nonlinearity (i.e. higherS/S and lowerm/m) and/or lower numbers of the nonlinear terms.
Before describing the transformations in a greater detail, the usual vector notation is introduced and related to the formal sum notation from Sect. 2. The vector notation is used when ODE systems are considered in matrix form, while the formal sum notation is used when ODE systems are considered component-wise.
Notation Let |S| = S, |C| = C and |R| = R, and suppose S, C and R are each given a fixed ordering with indices being n = 1, 2, . . . , S, i = 1, 2, . . . , C, and l = 1, 2, . . . , R, respectively, i.e. one can identify the ordered components of formal sums with components of Euclidean vectors. Let also the indices s n be denoted by n, n = 1, 2, . . . , S, for brevity. Then, the kinetic equations under the mass action kinetics in the formal sum notation are given by (1). In this section, we start with equations which have more general polynomial, and not necessarily kinetic, functions on the RHS, i.e. the ODE system is written in the formal sum notation as (8), while in the usual vector notation by where P(x; k) ∈ P m (R S ; R S ), x ∈ R S ≥ , and k ∈ R R .

Affine transformation
Definition 3.4 Consider applying an arbitrary nonsingular matrix A ∈ R S×S on Eq. (9), resulting in: wherex = Ax, andk is a vector of new rate constants obtained from k by rewriting the polynomial on the RHS of (10) into the mass action form. Then Ψ A : P m (R S ; R S ) → Definition 3.5 Consider substitutingx = x + T in Eq. (9), where T ∈ R S , which results in: wherek is a vector of the new rate constants obtained from k by rewriting the polynomial on the RHS of (11) into the mass action form.
A composition of a translation and a centroaffine transformation, i.e. an affine transformation, may be used as a possible kinetic transformation (see Definition 3.2). Let us note that condition (ii) in Definition 3.2 is necessarily satisfied for all affine transformation, i.e. affine transformations preserve the topology of the phase space, as well as the polynomial degree of the functions being mapped [25]. For these reasons, affine transformations are preferred over the alternative nonlinear transformations, discussed in the next two sections. However, affine transformations do not necessarily satisfy condition (i) in Definition 3.2, so that they are generally nonkinetic transformations. However, despite being generally nonkinetic, affine transformations map sets k into new setsk [see Eqs. (10), (11)], so that they may be used for satisfying a given set of constraints imposed on the output function (see Definition 3.3). This motivates the following definition.
is a kinetic function, and a given set of constraints Φ ≡ {φ j (k) ≥ 0 : j = 1, 2, . . . , J } is satisfied, then it is said that P(x; k) and the corresponding Eq. (8) are affinely kinetic, given the constraints. Otherwise, they are said to be affinely nonkinetic, given the constraints.
If the set of constraints in Definition 3.3 is empty, affinely nonkinetic functions are called essentially nonkinetic, while those that are affinely kinetic are called removably nonkinetic. Such labels emphasize that, if a function is essentially nonkinetic, a kinetic function that is globally topologically equivalent cannot be obtained, while if a function is removably nonkinetic, a globally topologically equivalent kinetic function can be obtained. Explicit sufficient conditions for a polynomial function P(x; k) to be affinely kinetic, or nonkinetic, are generally difficult to obtain. Even in the simpler case P(x; k) ∈ P 2 (R 2 ; R 2 ), such conditions are complicated, and cannot be easily generalized for higher-dimensional systems and/or systems with higher degree of nonlinearity [25]. In [19], based on the polar and spectral decomposition theorems, it has been argued that if no orthogonal transformation is kinetic, then no centroaffine transformation is kinetic. The result is reproduced in this paper using the more concise singular value decomposition theorem, and is generalized to the case when the set of constraints is nonempty. Loosely speaking, the theorem states that "orthogonally nonkinetic" functions are affinely nonkinetic as well, given certain constraints.
Proof By the singular value decomposition theorem, nonsingular matrices A ∈ R S×S can be written as A = Q 1 ΛQ 2 , where Q 1 , Q 2 ∈ R S×S are orthogonal, and Λ ∈ R S×S diagonal and positive definite. Cross-negative terms are invariant under transformation Ψ Λ for all Λ [19]. If Φ from Definition 3.3 is such that functions sign(φ j (k)), j = 1, 2, . . . , J , are invariant under all positive definite diagonal matrices Λ ∈ R S×S , the statement of the theorem follows.

X-factorable transformation
Then called an x-factorable transformation. If X is diagonal and its nonzero elements are where S ⊂ S, S = ∅, then the transformation is denoted Ψ X S , and is said to be x S -factorable.
Functions P(x; k) and (Ψ X P)(x; k) are not necessarily topologically equivalent due to two overlapping artefacts that Ψ X can produce, so that Ψ X is generally a nonkinetic transformation. Firstly, the fixed points of the former system can change the type and/or stability under Ψ X , and, secondly, the latter system has an additional finite number of boundary fixed points. The following theorem specifies the details of the artefacts for two-dimensional systems.

Theorem 3.3
Let us consider the ODE system (9) in two dimensions with RHS P(x; k) = (P 1 (x; k), P 2 (x; k)) . The following statements are true for all the fixed points x * of the two-dimensional system (9) in R 2 > under Ψ X S , S ⊆ S, S = ∅: (i) All the saddle fixed points are unconditionally invariant, i.e. saddle points of (9) correspond to saddle points of (12). (ii) A sufficient condition for stability of a fixed point x * to be invariant is: (iii) A sufficient condition for the type of a fixed point x * to be invariant is: Assume that the ODE system (9) does not have fixed points on the axes of the phase space. Nevertheless, the two-dimensional system (12) can have additional fixed points on the axes of the phase space, called boundary fixed points, denoted x * b ∈ R 2 ≥ . The boundary fixed points can be either nodes or saddles, and the following statements are true: , and the corresponding eigenvectors along the phase space axes.

a boundary fixed point is a node if and only if
with the node being stable if P i (x * b ; k) < 0, and unstable if P i (x * b ; k) > 0, i, j = 1, 2, i = j. Otherwise, the fixed point is a saddle.
Proof Without loss of generality, we consider two forms of the system (12) with S = 2: where p ∈ {0, 1}, so that system (13)- (14) is x-factorable for p = 1, but only x 1factorable for p = 0. The results derived for an x 1 -factorable system hold when the system is x 2 -factorable, if the indices are swapped. By writing P i (x; k) = P i , i = 1, 2, the Jacobian of (13)- (14), J X , is for p ∈ {0, 1} given by First, consider how fixed points of P(x; k) are affected by transformation Ψ X S . Denoting the Jacobian of two-dimensional system (9) by J , and assuming the fixed points are not on the axes of the phase space (i.e. x * ∈ R 2 > ), the Jacobians evaluated at x * are given by: Comparing the trace, determinant and discriminant of J (x * ) and J X (x * ), we deduce (i)-(iii).
To prove (iv)-(v), we evaluate J X at the boundary fixed points of the form If p = 1, then one of the boundary fixed points is x * b = 0, and the Jacobian becomes a diagonal matrix, so that condition (iv) holds. If x * b,2 = 0, then P 2 (0, x * b,2 ; k) = 0 in (15), and comparing the trace, determinant and discriminant of J (x * ) and J X (x * b ), we deduce (v). Theorem 3.3 can be used to find conditions that an x-factorable transformation given by Ψ X : P m (R 2 ; R 2 ) → P m+1 (R 2 ; R 2 ) is a kinetic transformation. While conditions (ii)-(iii) in Theorem 3.3 may be violated when Ψ X is used, so that Ψ X is a nonkinetic transformation, a composition of an affine transformation and an x-factorable transformation, i.e. Ψ X ,A,T = Ψ X • Ψ A • Ψ T , may be kinetic. Furthermore, such a composite transformation may also be used to control the boundary fixed points introduced by Ψ X . Finding an appropriate A and T to ensure the topological equivalence near the fixed points typically means that the region of interest in the phase space has to be positioned at a sufficient distance from the axes. However, since the introduced boundary fixed points may be saddles, this implies that the phase curves may be significantly globally changed, regardless of how far away from the axes they are. The most desirable outcome of controlling the boundary fixed points is to eliminate them, or shift them outside of the nonnegative orthant. The former can be attempted by ensuring that the nullclines of the original ODE system (9) do not intersect the axes of the phase space, while the latter by using the Routh-Hurwitz theorem [26].
An alternative transformation, which is always kinetic, that also does not change the dimension of an ODE system is the time-parametrization transformation [27]. However, while Ψ X increases the polynomial degree by one, and introduces only a finite number of boundary fixed points which are given as solutions of suitable polynomials, the time-parameterization transformation generally increases the nonlinearity degree more than Ψ X , and introduces infinitely many boundary fixed points.

The quasi-steady state transformation
The quasi-steady state assumption (QSSA) is a popular constructive method for reducing the dimension of ODE systems by assuming that, at a given time-scale, some of the species reach a quasi-steady state, so that they can be described by algebraic, rather than differential equations. The QSSA is based on Tikhonov's theorem [28,29] that specifies conditions ensuring that the solutions of the reduced system are asymptotically equivalent to the solutions of the original system. The original system is referred to as the total system, and it consists of the reduced subsystem, referred to as the degenerate system, and the remaining subsystem, called the adjoined system, so that the QSSA consists of replacing the total system with the degenerate one, by eliminating the adjoined system. Korzukhin's theorem [28,29] is an existence result ensuring that, given any polynomial degenerate system, there exists a corresponding total system that is kinetic.
Thus, Tikhonov's theorem can be seen as a constructive direct asymptotic dimension reduction procedure, while Korzukhin's theorem as an inverse asymptotic dimension expansion existence result. Korzukhin's theorem has an important implication that an application of the QSSA can result in a degenerate system that is structurally different than the corresponding total system. In this paper, the QSSA is assumed to necessarily be compatible with Tikhonov's theorem. If this is not the case, then it has been demonstrated in [27,30] that application of a QSSA can create dynamical artefacts, i.e. it can result in degenerate systems, not only structurally different, but also dynamically different from the total systems. The artefacts commonly occur due to the asymptotic parameters in Tikhonov's theorem not being sufficiently small. For example, it has been shown that exotic phenomena such as multistability and oscillations can exist in a degenerate system, while not existing in the corresponding total system [27,30].
Using Korzukhin's and Tikhonov's theorems, a family of kinetic total systems for an arbitrary nonkinetic polynomial degenerate system can be constructed, as is now shown. For simplicity, we denote x c = s∈S x c s s , for any x ∈ R S and c ∈ N S . Definition 3.8 Consider Eq. (8), and assume that the reaction set is partitioned, R = R 1 ∪ R 2 , R 1 ∩ R 2 = ∅, so that (8), together with the initial conditions, can be written as where x ∈ R S , α sr ∈ N S , β sr ∈ N S , α sr = β sr , a sr ∈ R ≥ and b sr ∈ R ≥ for all s ∈ S and r ∈ R. Assume further that the species set is partitioned, S = S 1 ∪ S 2 , S 1 ∩ S 2 = ∅, so that equations for species s ∈ S 1 are kinetic, while those for species s ∈ S 2 are nonkinetic. Consider the following total system, consisting of a degenerate system given by which satisfies the initial condition (17) with x 0 s > 0 for s ∈ S 2 , and an adjoined system given by where μ > 0, ω s > 0 are parameters, and p(x) is a polynomial function of x , mapping the RHS of differential equations in system (16)-(17), denoted P(x; k), to the RHS of differential equations of system (18)-(21), denoted (Ψ QSSA P)({x, y};k), with the constraint that x s > 0 for s ∈ S 2 , is called a quasisteady state transformation. Here,m ≤ m + m 0 + 2, andk is a vector of the new rate constants obtained from k by rewriting the polynomial (Ψ QSSA P)({x, y};k) into the mass action form. (16)-(17) and (18)-(21), corresponding to P(x; k) and (Ψ QSSA P)({x, y};k), are asymptotically equivalent in the limit μ → 0, and (Ψ QSSA P)({x, y};k) is a kinetic function.

Theorem 3.4 Solutions of systems
Proof Fixed points of system (20) satisfy y * s = ω s (x s p s (x)) −1 . The fixed points are isolated, and, since (from Definition 3.8) x s > 0 and p s (x) > 0, ∀x ∈ R S ≥ , ∀s ∈ S 2 , it follows that the fixed points are globally stable. Thus, the conditions of Tikhonov's theorem [28] are satisfied by the total system (18)- (21). Then, by applying the theorem, i.e. substituting y * s into (19), one recovers the corresponding degenerate system given by (16)- (17). Finally, the total system (18)-(21) is kinetic, as can be verified by using Definition 2.7.

Corollary 3.1
The quasi-steady state transformation Ψ QSSA , defined in Definition 3.8, is a kinetic transformation in the limit μ → 0.
An alternative transformation, for which condition (i) in Definition 3.2 is also always satisfied, and that also expands the dimension of an ODE system, is an incomplete Carleman embedding [2,31]. However, condition (ii) in Definition 3.2 is satisfied for the incomplete Carleman embedding only provided initial conditions for the adjoined system are appropriately constrained, and, furthermore, the transformation generally results in an adjoined system with a higher nonlinearity degree when compared to Ψ QSSA . In fact, Theorem 3.4 can be seen as an asymptotic alternative to the incomplete Carleman embedding, i.e. instead of requiring adjoined variables to satisfy y s (t) = ω s x −1 s (t) p −1 s (x(t)), one requires them to satify lim μ→0 y s (t) = ω s x −1 s (t) p −1 s (x(t)). The theorem can also be seen as an extension of using the QSSA to represent reactions of more than two molecules as a limiting case of bimolecular reactions [32] to the case of using the QSSA to represent cross-negative terms as a limiting case of kinetic ones.

Construction of reaction systems undergoing a supercritical homoclinic bifurcation
In this section, a brief review of a general bifurcation theory, and a more specific homoclinic bifurcation, is presented. This is followed by applying the framework developed in Sect. 3 to construct specific reaction systems displaying the homoclinic bifurcation.
Variations of parameters in a parameter dependent ODE system may change topology of the phase space, e.g. a change may occur in the number of invariant sets or their stability, shape of their region of attraction or their relative position. At values of the bifurcation parameters at which the system becomes topologically nonequivalent it is said that a bifurcation occurs, and the bifurcation is characterized by two sets of conditions: bifurcation conditions defining the type of bifurcation, and genericity conditions ensuring that the system is generic, i.e. can be simplified near the bifurcation to a normal form [33]. If it is sufficient to analyse a small neighbourhood of an invariant set to detect a bifurcation, the bifurcation is said to be local. Otherwise, it is called global, and the analysis becomes more challenging. Bifurcations are common in kinetic equations, where, in the case of the mass action kinetics, the rate constants play the role of bifurcation parameters [4,5,7,34,35]. In this paper, focus is placed on a global bifurcation giving rise to stable oscillations, called a supercritical homoclinic bifurcation [4,33,36]. Definition 4.1 Suppose x * is a fixed point of system (9). An orbit γ * starting at a point x is called homoclinic to the fixed point x * if its α-limiting and ω-limiting sets are both equal to x * .
Put more simply, a homoclinic orbit connects a fixed point to itself. An example of a homoclinic orbit to a saddle fixed point can be seen in Fig. 1b on page 20, where the homoclinic orbit is shown as the purple loop, while the saddle as the blue dot at the origin.
If a homoclinic orbit to a hyperbolic fixed point is present in an ODE system, then the system is structurally unstable, i.e. small perturbations to the equations can destroy the homoclinic orbit and change the structure of the phase space, so that a bifurcation can occur. For two-dimensional systems, the bifurcation and genericity conditions are completely specified by the Andronov-Leontovich theorem [33] given in Appendix 2. In summary, the theorem demands that the sum of the eigenvalues corresponding to the saddle at the bifurcation point, called the saddle quantity, must be nonzero (nondegeneracy condition), and that the so-called Melnikov integral at the bifurcation point evaluated along the homoclinic orbit must be nonzero (transversality condition).

The inverse problem formulation
Construction of reaction systems with prescribed properties is an inverse problem which we will solve by applying kinetic transformations described in Sect. 3. Our goal is to find a reaction system with the mass action kinetics (see Definition 2.6)

Construction of a polynomial function P(x; k):
Find an ODE system (8) which satisfies the assumptions of Andronov-Leontovich theorem in Appendix 2.

Construction of a kinetic function K(x;k):
Find a transformation so that the following conditions are satisfied: (i) The transformation is kinetic (see Definition 3.2), mapping the polynomial function P(x; k) into a kinetic function K(x;k) ≡ ( P)(x;k). (ii) The set of constraints (see Definition 3.3) ensuring that the homoclinic orbit is in R 2 ≥ are satisfied for K(x;k).
To determine the choice of , if possible, use Theorem 3.1 to deduce that P(x; k) is affinely nonkinetic (see Definition 3.6), given the constraints.

Construction of a reaction network:
Use Definition 3.1 to construct the canonical reaction network R K −1 .

Algorithm 1 Three steps of the solution to the inverse problem of finding reaction systems undergoing a supercritical homoclinic bifurcation.
such that the kinetic equations satisfy assumptions of Andronov-Leontovich theorem in Appendix 2, i.e. they must contain a homoclinic orbit defined on a two-dimensional manifold in the nonnegative orthant, and undergo the homoclinic bifurcation in a generic way. The output of this inverse problem will be a canonical reaction network which corresponds to the constructed ODE system. Thus the inverse problem is solved in three steps given in Algorithm 1. The first step is solved by using results of Sandstede [36] which leads to a set of polynomial functions satisfying the first three conditions of Andronov-Leontovich theorem in Appendix 2. An additional transformation is then performed ensuring that the final condition of Andronov-Leontovich theorem is satisfied. In this paper, nonlinear kinetic transformations are applied on the resulting polynomial function (using Step (2), case (b), in Algorithm 1).

Step (1): construction of polynomial function P(x; k)
Definition 4.2 A version of a plane algebraic curve called Tschirnhausen cubic (also known as Catalan's trisectrix, and l'Hospital's cubic) [37] given by: is referred to as the alpha curve. The part of the curve with x 2 < 0 is called the alpha loop, while the part with x 2 > 0 is called the alpha branches, with the branch for x 1 < 0 being the negative alpha branch, and for x 1 > 0 being the positive alpha branch. Solutions x 1 of Eq. (22) are denoted x ± 1 = ±x 2 The α curve is shown in Fig. 1a, with x − , and the positive and negative alpha branches in the first and second quadrant, respectively. As is done in [36], the alpha curve is mapped into a system of polynomial ODEs.

Lemma 4.1
The two-dimensional polynomial ODE system contains the alpha curve (22) as phase plane orbits, with the alpha loop as a homoclinic orbit to the fixed point x * = 0, provided a 2 < 1. If a ∈ (−1, 0), the alpha loop is stable from the inside, and system (23)-(24) has three fixed points: a saddle at the origin, an unstable spiral inside the alpha loop, and a stable node on the positive alpha branch.
A representative phase plane diagram of system (23)-(24) is shown in Fig. 1b. Note that a part of the positive alpha branch x + 1 is a heteroclinic orbit connecting the saddle and the node. The distance between the saddle and the node is given by d(a) = a −3 (a 2 − 1) √ 1 + a 2 , so that lim a→0 d(a) = +∞ and lim a→−1 d(a) = 0, i.e. increasing a increases length of the heteroclinic orbit by sliding the node along x + 1 . System (23)-(24) satisfies the first three conditions of Andronov-Leontovich theorem in Appendix 2. In order to satisfy the final condition, a set of perturbations must be found that destroy the alpha loop in a generic way, and this is ensured by the Melnikov condition (43). The bifurcation parameter controlling the existence of the alpha loop is denoted as α ∈ R. Note that P(x; a) perturbed by a function of the form α∇ H (x 1 , x 2 ) = α(−2x 1 , 2x 2 + 3x 2 2 ) satisfies the Melnikov condition [36], but P(x; a) + α∇ H (x 1 , x 2 ) has three terms dependent on α. In the following theorem, a simpler set of perturbations is found, introducing only one α dependent term in system (23)-(24). (α f (x 1 ), 0) T , where α ∈ R, is added to the RHS of system (23)- (24), and if f (x 1 ) is an odd function, f (−x 1 ) = − f (x 1 ), and f (x 1 ) = 0, ∀x 1 ∈ [−2 √ 3/9, 0), then the perturbed system undergoes a supercritical homoclinic bifurcation in a generic way as α is varied in the neighbourhood of zero.

Theorem 4.1 If a perturbation of the form
Proof Consider the perturbed version of system (23)-(24): Melnikov integral (43) for system (26)-(27) is given by Using (27), we have P 2 (x 1 , x 2 ; a)dt = dx 2 . Thus we can express M(0) in terms of x 2 as follow: where t + (x 2 ) (resp. t − (x 2 )) is the dependence of time t on x 2 along the positive (resp. negative) alpha branch for the trajectory which is at point (x 1 , x 2 ) = (0, −1) at time t = 0 (for α = 0). Since f is odd and ϕ(t ± ) > 0, we deduce For further simplicity of (26)- (27), function f (x 1 ) is set to f (x 1 ) = x 1 in the rest of this paper.

Step (2): construction of kinetic function K(x;k)
The RHS of system (26)- (27), P(x; a, α), is a kinetic function. However, the alpha loop, which is the region of interest, is not located in the nonnegative orthant. In order to position the loop into the positive orthant, we will apply affine transformations. First, we show that system (26)- (27) with the homoclinic orbit in nonnegative orthant is nonkinetic under all translation transformations for a ∈ (−1, 0), α ∈ R. P(x; a, α), given by the RHS of (26)- (27), is nonkinetic under all translation transformations Ψ T , for a ∈ (−1, 0) and α ∈ R, given the condition that the homoclinic orbit is mapped into R 2 > .
One can also readily prove that P(x; a, α) is nonkinetic under all affine transformations for |a| 1, and |α| 1. Thus, in the next two sections, we follow Step (2), case (b), in Algorithm 1, applying transformations Ψ X and Ψ QSSA on the kinetic function (Ψ T P)(x;k) given by (28). We require the following conditions to be satisfied: with the set of constraints Φ ensuring that the homoclinic orbit is in R 2 > . The reason for requiring |α| 1 is that then the following results are more readily derived, and the condition is sufficient for studying system (28) near the bifurcation point α = 0. A representative phase plane diagram corresponding to the ODE system with RHS (28) is shown in Fig. 2a, with fixed points, the alpha curve, and the vector field denoted as in Fig. 1b, and with the red segments on the axes corresponding to the phase plane regions where the cross-negative effect exists (see Definition 2.9).

X-factorable transformation
Let us apply the x-factorable transformation Ψ X on system (28). Letting Ψ X ,T ≡ Ψ X • Ψ T , the resulting kinetic function K X ,T (x;k) ≡ (Ψ X ,T P)(x;k) is given by , T 2 ∈ max(1, −aT 1 ), Proof Let us assume α = 0. Phase plane diagrams of a ODE system with RHS (28); b ODE system with RHS (32); and c ODE system with RHS (37). The stable node, the saddle, and the unstable spiral are represented as the green, blue and red dot, respectively. The alpha curve is shown in purple, and the vector field as gray arrows. On each plot it is indicated which kinetic transformation is applied to system (26)- (27). Red segments of the phase plane axes in (a) denote the regions where the cross-negative effect exists. In b and c, boundary fixed points are represented as the black dots, purple curves with a coarser dashing represent the saddle manifolds that asymptotically approach an axis, while with a finer dashing those that are outside of R 2 ≥ . d Phase plane diagram of a system for which x-factorable transformation significantly globally influences the phase curves in such a way that a pure translation cannot resolve the artefacts. For more details, see the text. The parameters are fixed to a = −0.8, α = 0, T 1 = 2, T 2 = 2 (Color figure online) From statement (i) of Theorem 3.3 it follows that the saddle fixed point of (28) is preserved under Ψ X . Denoting the node and spiral fixed points of (28) byx * nd andx * sp , respectively, one finds that the Jacobian is given by: Conditions (ii) and (iii) of Theorem 3.3 are both satisfied for the node, but only condition (ii) is satisfied for the spiral. The condition for the spiral to remain invariant is obtained by demanding disc(J X ,T |x =x * sp ) < 0, where J X ,T is the Jacobian of (32), and, taking into consideration (31), this leads to Boundary fixed points are given by (0, 0), (T 1 + a −1 T 2 , 0) and (0, 1/2(1 ± 1 + 4a −1 T 1 ) + T 2 ). The second fixed point can be removed from the nonnegative quadrant by demanding T 2 > −aT 1 , while the pair of the last ones always has nonzero imaginary part due to (31). Statement (iv) of Theorem 3.3 implies that the eigenvalues at the origin are given by λ 1 = k 1 0 = 3/2(T 2 − 2/3)(aT 1 + T 2 ) > 0 and λ 2 = k 2 0 = −T 1 + aT 2 (T 2 − 1) < 0, so origin is a saddle fixed point. As α can be chosen arbitrarily close to zero, and as K X ,T (x;k) is a continuous function of α, the theorem holds for sufficiently small |α| = 0, as well.
A representative phase plane diagram corresponding to the ODE system with RHS (32) is shown in Fig. 2b, where the saddle fixed point at the origin is shown as the black dot. It can be seen that one of the stable manifolds of the nonboundary saddle, represented as a dashed purple curve, approaches x 2 -axis asymptotically, instead of crossing it as in Fig. 2a.
The homoclinic orbit of the ODE system with RHS (32) is positioned below the node in the phase plane. Suppose the relative position of the stable sets is reversed by, say, applying an improper orthogonal matrix with the angle fixed to 3π/2, Ψ Q 3π/2− , with a representative phase plane shown in Fig. 2d. In this case, one can straightforwardly show that the boundary fixed point given by T 1 + 1/2(1 + 1 − 4a −1 T 2 ), 0 , shown as the black dot in Fig. 2d, cannot be removed from R 2 ≥ , and is always a saddle. The same conclusions are true for other similar configurations of the stable sets obtained by rotations only. This demonstrates that x-factorable transformation can produce boundary artefacts that have a significant global influence on the phase curves, that cannot be eliminated by simply translating a region of interest sufficiently far away from the axes. In order to eliminate the particular boundary artefact, a shear transformation may be applied. Consider applying Ψ X ,T ,Q 3π/2− , and T 1 = T 2 ≡ T ∈ R, for simplicity, leading to K n,X ,T ,Q 3π/2− ,S 2 (x;k) =x n (k n 0 +k n 1x 1 +k n 2x 2 +k n 11 (x 1 ) 2 +k n 12x 1x2 +k n 22 (x 2 ) 2 ), n = 1, 2, (37) where the coefficientsk =k(a, α, T ) are given by: , together with a ∈ (−1, 0), |α| 1 and Φ = {T > 1}. (26)- (27) and (37) are topologically equivalent in the neighbourhood of the fixed points in R 2 > , with the homoclinic orbit in R 2 > , and a saddle at the origin and a saddle with coordinates (0, 1/2T a −1 (1 + 2a)) being the only boundary fixed points in R 2 ≥ , if:
Note that as α → − 1 2 , the only boundary fixed point in R 2 ≥ is a saddle at the origin, and it is connected via a heteroclinic orbit to the saddle in R 2 > . A representative phase plane diagram corresponding to the ODE system with RHS (37) is shown in Fig. 2c.
While systems (32) and (37) contain specific variations of the specific homoclinic orbit given by (22), they, nevertheless, indicate possible phase plane topologies of the kinetic equations containing homoclinic orbits of shapes similar to (22). With a fixed shape and orientation of a homoclinic loop which is similar to (22), three possible orientations of a corresponding saddle manifold in R 2 > are: it may extend in R 2 > without asymptotically approaching a phase plane axis, it may asymptotically approach an axis, or it may cross an axis at a fixed point. In Fig. 2b, a combination of the first and second orientation is displayed, while in Fig. 2c of the first and third orientation.
Note that a composition an x-factorable and a quasi-steady state transformation may be used to make (28) kinetic. For example, one may eliminate the cross-negative termk 1 2 in (28) by using the x 1 -factorable transformation Ψ X 1 , and the cross-negative termk 2 0 by using an appropriate Ψ QSSA . An example of a kinetic function obtained by a transformation of the form Ψ QSSA,X 1 ,T is given by withx 2 (0) > 0, ω > 0, and μ → 0. Note that the chosen Ψ QSSA,X 1 ,T does not introduce any additional fixed points when applied to system (28).
leads to autocatalytic reaction networks, has been defined in Sect. 3.2.2, and its properties when applied to two-dimensional systems have been derived in Theorem 3.3. The quasi-steady state transformation, that increases the dimension of the systems being transformed, but generally introduces a lower number of nonlinear terms, has been presented in Sect. 3.2.3, and justified using Tikhonov's and Korzukhin's theorems [28]. As the focus of the paper has been more placed on the construction of kinetic equations, and less on constructions of reaction networks, in Sect. 3.1 an analytical and algorithmic method for obtaining the so-called canonical networks has been presented [6]. The framework may be used for constructing lower-dimensional reaction systems displaying exotic phenomena, with Algorithm 1 exemplifying the construction process.
In the second part of the paper, the inverse problem framework has been applied to a case study of constructing bistable reaction systems undergoing a supercritical homoclinic bifurcation, with a parameter controlling the stable sets separation, with an overview of the procedure presented in Sect. 4.1. In Sect. 4.2, a polynomial ODE system having a homoclinic orbit in the phase plane has been constructed using the results from [36], and perturbed in such a way that the sufficient conditions for the existence of the homoclinic bifurcation are fulfilled, as required by Andronov-Leontovich Theorem [33]. In Sect. 4.3, the kinetic transformations from Sect. 3 have been used in order to map the polynomial system to a kinetic one with the homoclinic orbit in the positive quadrant. The topological phase space effects produced by the kinetic transformations on the constructed systems have been discussed. In Sect. 4.4 and Appendix 3, the canonical reaction networks induced by some of the kinetic equations have been presented.
In this paper, we have constructed chemical reaction networks inducing twodimensional cubic kinetic functions with the deterministic ODEs (kinetic equations) undergoing a supercritical homoclinic bifurcation. In a future publication, we will report our results on the stochastic analysis of the constructed systems, consisting of examining the quasi-stability of the limit cycle, and stochastic switching between the stable sets, as a function of the bifurcation parameter and the parameter controlling the stable set separation. A motivation for such a study is the fact that stochastic effects play an important role in systems biology due to the inherently small reactors [5,7,34,35], and one might even say that systems biology has initiated a renaissance of the stochastic reaction kinetics [38].
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Appendix 1: Oscillations in a subset of two-dimensional second-order reaction systems with cross-negative terms Theorem 5.1 Consider a two-dimensional ODE system with a quadratic polynomial RHS, P(x; k) ∈ P 2 (R 2 ; R 2 ), k ∈ R 12 : dx 1 dt = P 1 (x; k) = k 1 0 + k 1 1 x 1 + k 1 2 x 2 + k 1 11 x 2 1 + k 1 12 x 1 x 2 + k 1 22 x 2 2 , dx 2 dt = P 2 (x; k) = k 2 0 + k 2 1 x 1 + k 2 2 x 2 + k 2 11 x 2 1 + k 2 12 x 1 x 2 + k 2 22 If k 1 11 ≤ 0 and k 2 22 ≤ 0, and if system (42) is nonnegative, then the system has no limit cycles.
Then, for all sufficiently small |α|, there exists a neighbourhood of the saddle fixed point and the homoclinic orbit such that a unique limit cycle bifurcates from γ * . If σ 0 < 0, the homoclinic bifurcation is supercritical, giving rise to a unique stable limit cycle, while if σ 0 > 0, the homoclinic bifurcation is subcritical, giving rise to a unique unstable limit cycle.