CHEMICAL REACTION SYSTEMS WITH A HOMOCLINIC BIFURCATION: AN INVERSE PROBLEM

: A framework for inverse problems is offered to build reaction systems with specified features. The framework includes the definition and analysis of kinetic transformations, which enable the mapping of any polynomial ordinary differential equation to the one that may be represented as a reaction network. The framework is applied to the design of certain bistable reaction systems in two and three dimensions that experience a supercritical homoclinic bifurcation, and the phase spaces' topology is examined.


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 analytic ODEs, thus being important for pure mathematics [16,18].They can display not only the chemically regular phenomenon of having a globally stable fixed point, but also the chemically exotic phenomena (multistability, limit cycles, 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 [23].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 [26,27].The construction of reaction networks displaying prescribed properties may be seen as an inverse problem in formal reaction kinetics [4], where, given a set of properties, a set of compatible reaction networks is searched for.Such constructions are useful in application areas such systems biology (as caricature models), synthetic biology (as blueprints), and numerical analysis (as test problems) [25,30].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 [30].A crucial property of the kinetic equations is a lack of so-called cross-negative terms [2], 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 existence of cross-negative terms, together with a requirement that the dependent variables are always finite, imply that not every nonnegative polynomial ODE system is kinetic, and, thus, further constrain the possible dynamics, playing an important role in the construction of reaction systems, chemical chaos, and pattern formation via Turing instabilities [3,5].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 limit 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 constrained two-dimensional kinetic ODE systems.Motivated by this, this paper consists of two main results: firstly, building upon the framework from [4,2], an inverse problem framework suitable for constructing the reaction systems is presented in Section 3, with the focus on the so-called kinetic transformations, allowing one to map a nonkinetic into a kinetic system.Secondly, in Section 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 Section 2. A summary of the paper is presented in Section 5.

Notation and definitions
The notation and definitions in this paper are inspired by [1,8,2].Definition 2.1.Let R be the space of real numbers, R≥ the space of nonnegative real numbers, R> the space of positive real numbers and N = {0, 1, 2, 3, . . .} the set of natural numbers.Given a finite set S, with cardinality |S| = S, the real space of formal sums c = P s∈S css is denoted by R S if cs ∈ R for all s ∈ S. It is denoted by R S ≥ if cs ∈ R≥ for all s ∈ S; by R S > if cs ∈ R> for all s ∈ S; and by N S if cs ∈ N for all s ∈ S; where the number cs is called the s-component of c for s ∈ S. Support of c ∈ R S is defined as supp(c) = {s ∈ S : cs 6= 0}.Complement of a set M ⊂ S is denoted by Mc , and given by Mc = S \ M. 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 equation (3.3).

Reaction networks and reaction systems
Elements r = c → c 0 are called reactions of the network, and it is said that c reacts to c 0 , with c being called the reactant complex, and c 0 the product complex.The order of reaction r is given by or = P s∈S cs < ∞ for r = c → c 0 ∈ 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 [2], in this paper we consider reactions with or ≤ 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 0 ) ∈ R, positivity κr(x) > 0 is satisfied if and only if supp(c) ⊂ supp(x).

Research Article
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) : R S ≥ → R S , is given by K(x; R) = P r∈R κr(x)(c 0 − c) where r = c → c 0 .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 (2.1) are required to be finite, i.e. xs < ∞, for s ∈ S, and for t ≥ 0, except possibly for initial conditions located on a finite number of (S − z)dimensional subspaces of R S ≥, z ≥ 1.

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.
Definition 2.7.Let f : R S ≥ → R S be given by fs(x) = P r∈R fsr(x), where fsr(x) ∈ R, for x ∈ R S ≥, s ∈ S and r ∈ R. If ∃s ∈ S, ∃r ∈ R and ∃x ∈ R S ≥ such that s ∈ suppc (x) and fsr(x) < 0, then fsr(x) is called a cross-negative term, and function f(x) and ODE system dx/dt = f(x) are said to be nonkinetic.

Research Article
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.Note that a system {R, k}, corresponding to N (x; k) in Definition 2.8, has a well-defined reaction network R (for r = c → c 0 , r ∈ R, we restrict c, c0 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.2).Note also that Pm(R S ≥; R S ) = P K m(R S ≥; R S )∪ P N m(R S ≥; R S ), with P K m(R S ≥; R S ) ∩ P N m(R S ≥; R S ) = ∅.

Properties of kinetic functions
From Definition 2.3 it follows that a kinetic function K(x; R) has a structural property: crossnegative 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.Definition 2.9.Let f : R S ≥ → R S be given by fs(x) = P r∈R fsr(x), where fsr(x) ∈ R, for x ∈ R S ≥, s ∈ S and r ∈ R. If ∀s ∈ S, ∀x ∈ R S ≥, s ∈ suppc (x) ⇒ fs(x) ≥ 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 ∈ suppc (x) and fs(x) < 0.
Note that the absence of cross-negative terms implies nonnegativity, but the converse is not necessarily true [2,7], i.e. an ODE system may have cross-negative terms, without having a crossnegative effect, as we will show in Example 2.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 [3,5].In the context of oscillations, as a generalization of the result in [20], one can prove that in twodimensional reaction systems with mass action form and with at most bimolecular reactions, the nonexistence of a cross-negative effect in the ODEs is a sufficient condition for nonexistence of limit cycles (see A). ) induces a reaction system only in case (i).In particular, nonnegative ODE system (2.5)-(2.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 (2.5)-(2.6) is nonnegative for all forward times).
Let f : R S ≥ → R S be given by fs(x) = P r∈R fsr(x), where fsr(x) ∈ R, for x ∈ R S ≥, s ∈ S and r ∈ R.  [10] (see also Section 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 (2.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 [2,4].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 Section 3.1, while the first step in Section 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.
3.1 The canonical reaction network Let us assume that we are able to construct an ODE system of the form (2.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 [3,4].While, for a fixed kinetics, a reaction network induces 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 [8,9,3].For example, in [9], 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 [3], 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.The procedure is used in this paper, so it is now defined.
Definition 3.1.Let κ : R S ≥ → R R ≥ be a kinetics.Consider the kinetic function given by K(x; R) = P r∈R drκr(x), where x ∈ R S ≥ and dr ∈ R S .Let us map K(x; R) to a reaction system {RK−1 , κK−1 } with complexes and kinetics given by: (i) Reactant complexes, cr, are assumed to be uniquely obtainable from κr(x) for r ∈ R, which is true in the case of the mass action kinetics.(ii) Reaction cr → c 0 rs is then constructed for each r ∈ R and s ∈ S, where new product complexes are given by c 0 rs = cr + sign(drs)s, with sign(•) being the sign function.(iii) The new kinetics is then defined as κK−1r 0(x) ≡ |drs|κr(x), for r ∈ R, s ∈ S, where r 0 ∈ RK−1 .The induced reaction system {RK−1 , κK−1 } is called the canonical reaction system, with RK−1 being the canonical reaction network.
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 straightforward way, while they generally do not contain minimal number of reactions.

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) [2].However, only kinetic functions induce reaction networks, as exemplified in Example 2.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.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 locally topologically equivalent to the input function in the neighbourhood of the corresponding fixed points, and the dynamics along any additional dimensions of the output function (corresponding to the additional species) must asymptotically tend to the corresponding fixed point.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.A subset of constraints is now defined.
There are two sets of kinetic transformations.The first, and the preferred, set of possible kinetic transformations are affine transformations, which are discussed in Section 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 cross-negative terms, and two such transformations are discussed in Sections 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.lower S/S ¯ and higher ¯m/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. higher S/S ¯ and lower ¯m/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 Section 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.A composition of a translation and a centroaffine transformation, ΨA,T = ΨA • ΨT , 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 [6].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 sets ¯k (see equations (3.4) and (3.5)), 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.
Definition 3.6.Let P(x; k) ∈ Pm(R S ; R S ).If it is not possible that simultaneously (ΨA • ΨT P)(¯x; ¯k) is a kinetic function, and that a given set of constraints Φ ≡ {φj ( ¯k) ≥ 0 : j = 1, 2, . . ., J} is satisfied, for all A ∈ R S×S and for all T ∈ R S , then it is said that P(x; k) and the corresponding equation (3.2) are affinely nonkinetic, given the constraints.Otherwise, they are said to be affinely kinetic, 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) ∈ P2(R 2 ; R 2 ), such conditions are complicated, and cannot be easily generalized for higher-dimensional systems and/or systems with higher degree of nonlinearily [6].In [5], 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.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 (3.3) in two dimensions with RHS P(x; k) = (P1(x; k),P2(x; k))>.The following statements are true for all the fixed points x * of the twodimensional system (3.3) in R 2 > under ΨXS0 , S 0 ⊆ S, S 0 6= ∅: (i) All the saddle fixed points are unconditionally invariant, i.e. saddle points of (3.3) correspond to saddle points of (3.6).(ii) (ii) A sufficient condition for stability of a fixed point x * to be invariant is: 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 6= 0, then P2(0, x * b,2 ; k) = 0 in (3.9), and comparing the trace, determinant and discriminant of J(x * ) and JX (x * b ), we deduce (v).Theorem 3.3 can be used to find conditions that an x-factorable transformation given by ΨX : Pm(R 2 ; R 2 ) → Pm+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 (3.3)do not intersect the axes of the phase space, while the latter by using the Routh-Hurwitz theorem [34].
An alternative transformation, which is always kinetic, that also does not change the dimension of an ODE system is the time-parametrization transformation [14].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 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 [11,12] 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 adjointed system, so that the QSSA consists of replacing the total system with the degenerate one, by eliminating the adjointed system.Korzukhin's theorem [11,12] 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 [14,15] 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 [14,15].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 = Q s∈S x cs s , for any x ∈ R S and c ∈ N S .Proof.Fixed points of system (3.14)satisfy y * s = ωs(xsps(x))−1 .The fixed points are isolated, and, since (from Definition 3.8) xs > 0 and ps(x) > 0, ∀x ∈ R S ≥, ∀s ∈ S2, it follows that the fixed points are globally stable.Thus, the conditions of Tikhonov's theorem [11] are satisfied by the total system (3.12)-(3.15).Then, by applying the theorem, i.e. substituting y * s into (3.13),one recovers the corresponding degenerate system given by (3.10)- (3.11).Finally, the total system (3.12)-(3.15) 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 [18,17].
However, condition (ii) in Definition 3.2 is satisfied for the incomplete Carleman embedding only provided initial conditions for the adjointed system are appropriately constrained, and, furthermore, the transformation generally results in an adjointed 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 adjointed variables to satisfy ys(t) = ωs x −1 s (t) p −1 s (x(t)), one requires them to satify limµ→0 ys(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 [19] 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 Section 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 [24].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 [25,26,27,28,29].In this paper, focus is placed on a global bifurcation giving rise to stable oscillations, called a supercritical homoclinic bifurcation [24, 21, 26].
Definition 4.1.Suppose x * is a fixed point of system (3.3).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 Figure 1(b) on page 16, 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 twodimensional systems, the bifurcation and genericity conditions are completely specified by the Andronov-Leontovich theorem [24] given in B. 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 Section 3. Our goal is to find a reaction system with the mass action kinetics (see Definition 2.6) such that the kinetic equations satisfy assumptions of Andronov-Leontovich theorem in B, i.e. they must contain a homoclinic orbit defined on a twodimensional 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 3) for a = −0.8, with the stable node, the saddle, and the unstable spiral represented as the green, blue and red dots, respectively.The alpha curve is shown in purple, while the vector field as gray arrows.curve, and the vector field denoted as in Figure 1(b), 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 (4.Note that a composition an x-factorable and a quasi-steady state transformation may be used to make (4.7) kinetic.For example, one may eliminate the cross-negative term ¯k 1 2 in (4.7) by using the x1factorable transformation ΨX1 , and the cross-negative term ¯k 2 0 by using an appropriate ΨQSSA.An example of a kinetic function obtained by a transformation of the form ΨQSSA,X1,T is given by with ¯x2(0) > 0, ω > 0, and µ → 0. Note that the chosen ΨQSSA,X1,T does not introduce any additional fixed points when applied to system (4.7).

Summary
In the first part of the paper, a framework for constructing reaction systems having prescribed properties has been formulated as an inverse problem and presented in Section 3, relying on definitions introduced in Section 2. As a part of the framework, in Section 3.2, kinetic transformations have been defined that enable one to map an arbitrary polynomial ODE system with a set of constraints, possibly containing the cross-negative terms, into a kinetic one.Augmented with the results from [16], such transformations can be applied to nonpolynomial systems as well.Systems for which no affine transformation is kinetic have been defined as affinely nonkinetic in Section 3.2.1, to emphasize the fact that significant changes to their solutions are required.X-factorable transformation [10], that does not change the dimension of the systems being transformed, but introduces a higher number of nonlinear terms and leads to autocatalytic reaction networks, has been defined in Section 3.2.2, and its properties when applied on 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 Section 3.2.3, and justified using Tikhonov's and Korzukhin's theorems [11].As the focus of the paper has been more placed on the construction of kinetic equations, and less on constructions of reaction networks, in Section 3.1 an analytical and algorithmic method for obtaining the so-called canonical networks has been presented [4].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 Section 4.1.In Section 4.2, a polynomial ODE system having a homoclinic orbit in the phase plane has been constructed using the results from [21], and perturbed in such a way that the sufficient conditions for the existence of the homoclinic bifurcation are fulfilled, as required by AndronovLeontovich Theorem [24].In Section 4.3, the kinetic transformations from Section 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 Section 4.4 and C, the canonical reaction networks induced by some of the kinetic equations have been presented.
In this paper, we have constructed chemical reaction networks inducing two-dimensional 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 [25, 27, 28, 29], and one might even say that systems biology has initiated a renaissance of the stochastic reaction kinetics [31].

A Oscillations in two-component bimolecular reaction systems with cross-negative terms
Theorem A.1: Consider a two-dimensional ODE system with a quadratic polynomial RHS, P(x; k) ∈ P2(R 2 ; R 2 ), k ∈ R 12: If k 1 11 ≤ 0 and k 2 22 ≤ 0, and if system (A.1) is nonnegative, then the system has no limit cycles. .If k 1 11 ≤ 0 and k 2 22 ≤ 0, and if system (A.1) is nonnegative, so that P1(0, x2) ≥ 0 and P2(x1, 0) ≥ 0, then D¯ ≥ 0. The only case when a limit cycle may exist is if D¯ = 0 for all x1, x2 > 0, and in [20] it is shown that no limit cycles exist in this case.
In [20], it was shown that the absence of cross-negative terms in P(x; k) ∈ P2(R 2 ; R 2 ), with k 1 11 ≤ 0 and k 2 22 ≤ 0, implies the absence of limit cycles, i.e. one requires the more restrictive condition P(x; k) ∈ P K 2 (R 2 ≥; R 2 ).Theorem A.1 shows that, in fact, absence of the cross-negative effect in

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 S 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, c0 ), denoted r = c → c 0 , with the following properties:

Definition 2 . 6 .
Kinetics κ is called the mass action kinetics if κr(x) = krx c , for r = (c → c 0 ) ∈ R, where kr > 0 is the rate constant of reaction r, and x c = Q s∈S x cs s , with 00 = 1.A reaction system with the mass action kinetics is denoted {R, k}, and the corresponding kinetic function is denoted K(x; k) ≡ K(x; R) = P r∈R kr(c 0 − c)x c , where k ∈ R R > .A review of the mass action kinetics can be found in [32].In this paper, most of the results are stated with kinetics fixed to the mass action kinetics.This is not restrictive, as an arbitrary analytic function can always be reduced to a polynomial one [16].Example 2.1.Consider the following reaction network (consisting of one reaction) under the mass action kinetics: so that S = {s1, s2}, C = {s1 + s2, 2s2}, R = {s1 + s2 → 2s2} and k = {k1}.Concentration x ∈ R S ≥ has two components.To simplify our notation, we write x1 = xs1 , x2 = xs2 , and K1(x; k) = Ks1 (x; R), K2(x; k) = Ks2 (x; R).Then the induced system of kinetic equations is given by

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.Then, the set of functions P(x; k) is denoted by Pm(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 ).

Definition 3 . 2 .
Let P(x; k) ∈ Pm(R S ; R S ), k ∈ R R, i.e.P(x; k) is a polynomial function.Consider the corresponing 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: (i) Ψ : Pm(R S ; R S ) → P K m¯ (R S¯ ≥; R S¯ ), ¯m ≥ m, S¯ ≥ S, maps the polynomial function P(x; k) into a kinetic function K(¯x; ¯k) ≡ (Ψ(P))(¯x; ¯k) for ¯x ∈ R S¯ ≥ and ¯k ∈ R R¯ > .(ii) (ii) Let x * be a fixed point of (3.2) that is mapped by Ψ to fixed point ¯x * ∈ R S¯ ≥ of the system of kinetic equations (2.1) with K(¯x; ¯k) on its RHS.Let also the eigenvalues of the Jacobian matrix of P(x; k), J(x * ; k), denoted by λn, n = 1, 2, . . ., S, be mapped to the eigenvalues of Jacobian of K(¯x; ¯k), JΨ(¯x * ; ¯k), which are denoted by λ¯ n, n = 1, 2, . . ., S.Then, for every such fixed point x * it must be true that sign(λn) = sign(λ¯ n), n = 1, 2, . . ., S, and, if there are any additional eigenvalues λ¯ n, n = S, S + 1, . . ., S¯, they must satisfy sign(λ¯ n) < 0.

Notation.Definition 3 . 5 .
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 sn 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 (2.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 (3.2), while in the usual vector notation by where P(x; k) ∈ Pm(R S ; R S ), x ∈ R S ≥, and k ∈ R R. 3.2.1 Affine transformation Definition 3.4.Consider applying an arbitrary nonsingular matrix A ∈ R S×S on equation (3.3), resulting in: where ¯x = Ax, and k¯ is a vector of new rate constants obtained from k by rewriting the polynomial on the RHS of (3.4) into the mass action form.Then ΨA : Pm(R S ; R S ) → Pm(R S ; R S ), mapping P(x; k) to (ΨAP)(¯x; k¯), is called a centroaffine transformation.If A is an orthogonal matrix, then ΨA is called an orthogonal transformation.Consider substituting ¯x = x + T in equation (3.3),where T ∈ R S , which results in:where k¯ is a vector of the new rate constants obtained from k by rewriting the polynomial on the RHS of (3.5) into the mass action form.Then ΨT : Pm(R S ; R S ) → Pm(R S ; R S ), mapping P(x; k) to (ΨT P)(¯x; k¯), is called a translation transformation.

Figure 1 :
Figure 1: (a) The alpha curve (4.1), with branch x − 1 plotted as the solid purple curve, and branch x + 1 as the dashed green curve.(b) Phase plane diagram of system (4.2)-(4.3)for a = −0.8, with the stable node, the saddle, and the unstable spiral represented as the green, blue and red dots, respectively.The alpha curve is shown in purple, while the vector field as gray arrows.

Theorem 4 . 3 . 4 . 3 . 2
7).Letting ΨX,T ≡ ΨX • ΨT , the resulting kinetic function KX,T (¯x; k¯) ≡ (ΨX,T P)(¯x; k¯) is given by Theorem 4.2.ODE systems with RHSs (4.7) and(4.11)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 being the only boundary fixed point in R 2 ≥, if:Proof.Let us assume α = 0.From statement (i) of Theorem 3.3 it follows that the saddle fixed point of (4.7) is preserved under ΨX .Denoting the node and spiral fixed points of (4.7) by ¯x * nd and ¯x * 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(JX,T |¯x=¯x * sp ) < 0, where JX,T is the Jacobian of (4.11), and, taking into consideration (4.10), this leads to Boundary fixed points are given by (0, 0), (T1+a −1T2, 0), (0, 1/2(1± √ 1 + 4a−1T1)+T2).The second fixed point can be removed from the nonnegative quadrant by demanding T2 > −aT1, while the pair of the last ones always has nonzero imaginary part due to (4.10).Statement (iv) of Theorem 3.3 implies that the eigenvalues at the origin are given by λ1 = k 1 0 = 3/2(T2 − 2/3)(aT1 + T2) > 0 and λ2 = k 2 0 = −T1 + aT2(T2 − 1) < 0, so origin is a saddle fixed point.As α can be chosen arbitrarily close to zero, and as KX,T (¯x; k¯) is a continuous function of α, the theorem holds for sufficiently small |α| 6= 0, as well.A representative phase plane diagram corresponding to the ODE system with RHS (4.11) is shown in Figure2(b), 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 x2-axis asymptotically, instead of crossing it as in Figure2(a).The homoclinic orbit of the ODE system with RHS (4.11) 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, ΨQ3π/2− , with a representative phase plane shown in Figure2(d).In this case, one can straightforwardly show that the boundary fixed point given by T1 + 1/2(1 + √ 1 − 4a−1T2), 0 , shown as the black dot in Figure2(d), 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 ,Q3π/2−,S2 = ΨX • ΨT • ΨQ3π/2− • ΨS2 on (4.5)-(4.6),where and T1 = T2 ≡ T ∈ R, for simplicity, leading to where the coefficients k¯ = k¯(a, α, T ) are given by: ODE systems with RHSs (4.5)-(4.6)and (4.16) 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; 18) Proof.Following the same procedure as in the proof of Theorem 4.2, and noting that the saddle, node and spiral fixed points are given by T , T , T − 2a −2 (1 − a 2 ), T + a −3 (1 − a 2 ) , and T + Research Article2/9(3+a 2 ), T −2/9a , respectively, while the five boundary fixed points are (0, 0), (T + 1/2 5T a± p T (8a−1(1 − a 2) + 9T a 2) , 0), (0, 1/2T a −1 (1+2a)), (0, a−1 2/3+T (1+a) ), one finds (4.18).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 (4.16) is shown in Figure2(c).While systems (4.11) and (4.16) contain specific variations of the specific homoclinic orbit given by (4.1), they, nevertheless, indicate possible phase plane topologies of the kinetic containing homoclinic orbits of shapes similar to (4.1).With a fixed shape and orientation of a homoclinic loop which is similar to (4.1), 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 Figure2(b), a combination of the first and second orientation is displayed, while in Figure 2(c) of the first and third orientation.The quasi-steady state transformation In Lemma 4.2, it was demonstrated that system (4.5)-(4.6)has at least one cross-negative term under translation transformations.It can be readily shown that system (4.5)-(4.6) in fact has minimally two cross-negative terms under translation transformations, ¯k 1 2 < 0 and ¯k 2 0 < 0, and this is the case when a ∈ (−1, 0), Φ = {T1 ∈ (2√ 3/9, −T2a), T2 > 1}.Let us apply a quasisteady state transformation ΨQSSA on system (4.7) that eliminates the two cross-negative terms, i.e. two new variables are introduced, ¯y1, y¯2 ∈ R 2 >, and we take p1(¯x) = p2(¯x) = 1, in Definition 3.8.Letting ΨQSSA,T ≡ ΨQSSA • ΨT , the resulting kinetic function KQSSA,T ({x, ¯ y¯}; ¯k, ω, µ) ≡ ΨQSSA,T P)({x, ¯ y¯}; ¯k, ω, µ) is given b with ¯xn(0) > 0, ωn > 0, n = 1, 2, and µ → 0. In (4.19), limxn→0 limµ→0 yn = +∞, n = 1, 2, and a geometrical implication is that, say, the saddle manifold crossing the x2-axis in Figure 2(a), instead asymptotically approaches the y1-axis.The asymptotic behavior of the saddle manifolds is achieved by the additional (boundary) fixed points in (4.11) displayed in Figure 2(b), and by additional phase space dimensions in (4.19).

4. 4
Step (3): construction of the canonical reaction network RK−1 Definition 3.1 can be used to map the kinetic functions (4.11), (4.16),(4.19),and (4.20) to the canonical reaction networks RK−1 .This is illustrated for (4.11) in this section, and for (4.16) and (4.20) in C. For clarity, both the induced kinetic equations and the induced canonical reaction networks are presented.Note that the reaction networks are assumed to be taking place in an open reactor, and are not necessarily purely chemical in nature.Nevertheless, the non-chemical processes present in kinetic equations are represented as quasi-chemical reactions.Such reactions take form of input/output of chemical species, as well as containing quasi-species that are time-independent on a relevant time scale, so that their constant concentration is absorbed into a quasi-kinetics, leading to conservation laws not necessarily holding[1].Writing x ≡ ¯x, the induced kinetic equations for (4.11) are given by while the induced canonical reaction where k1 = | ¯k 1 0 |, k2 = | ¯k 2 0 |, k3 = | ¯k 1 1 |, k4 = | ¯k 2 1 |, k5 = | ¯k 1 2 |, k6 = | ¯k 2 2 |, k7 = | ¯k 1 12|, k8 = | ¯k 1 22|, k9 = | ¯k 2 22|, with the coefficients k¯ given by (4.8), and the conditions given by (4.12).