Chemical Reaction Systems with a Homoclinic
Bifurcation: an Inverse Problem
Tomislav Plesa Mathematical Institute, University of Oxford, Radcliffe Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, United Kingdom; e-mails: firstname.lastname@example.org, email@example.com Tomáš Vejchodský Institute of Mathematics, Czech Academy of Sciences, Žitná 25, Praha 1, 115 67, Czech Republic, e-mail: firstname.lastname@example.org Radek Erban
Abstract: 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.
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 . 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 , 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 .
A crucial property of the kinetic equations is a lack of so-called cross-negative terms , 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 , for constant , where the term , although a polynomial of degree zero, nevertheless cannot be directly represented by a reaction, and results in .
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.
2 Notation and definitions
Let be the space of real numbers, the space of nonnegative real numbers, the space of positive real numbers and the set of natural numbers. Given a finite set , with cardinality , the real space of formal sums is denoted by if for all . It is denoted by if for all ; by if for all ; and by if for all ; where the number is called the -component of for . Support of is defined as . Complement of a set is denoted by , and given by .
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).
2.1 Reaction networks and reaction systems
A reaction network is a triple , where
is a finite set, with elements called the species of the network.
is a finite set, with elements , called the complexes of the network, such that . Components of are called the stoichiometric coefficients.
is a binary relation with elements , denoted , with the following properties:
such that or .
Elements are called reactions of the network, and it is said that reacts to , with being called the reactant complex, and the product complex. The order of reaction is given by for .
Note that as set implies sets and , reaction networks are denoted with , for brevity. Also, as it is unlikely that a reaction between more than three reactants occurs , in this paper we consider reactions with . To represent some of the non-chemical processes as quasireactions, the zero complex is introduced, denoted with , with the property that , where is the empty set.
Let be a reaction network and let be a continuous function which maps (called “species concentrations”) into (called “reaction rates”). Then is said to be a kinetics for provided that, for all and for all , positivity is satisfied if and only if .
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.
A reaction network augmented with a kinetics is called a reaction system, and is denoted .
Given a reaction system , the induced kinetic function, , is given by where . The induced system of kinetic equations, describing the time evolution of species concentrations , 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. , for , and for , except possibly for initial conditions located on a finite number of -dimensional subspaces of , .
Kinetics is called the mass action kinetics if , for , where is the rate constant of reaction , and , with . A reaction system with the mass action kinetics is denoted , and the corresponding kinetic function is denoted , where .
A review of the mass action kinetics can be found in . 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 .
Consider the following reaction network (consisting of one reaction) under the mass action kinetics:
so that , , and . Concentration has two components. To simplify our notation, we write , , and , . Then the induced system of kinetic equations is given by
2.2 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.
Let be given by , where , for , and If , and such that and , then is called a cross-negative term, and function and ODE system 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.
Let , , be a polynomial function with polynomial degree , . Then, the set of functions is denoted by . If is a kinetic function, it is denoted by , , and the set of such functions is denoted by . If is a nonkinetic function, it is denoted by , , and the set of such functions with domain is denoted by , while with domain by .
Note that a system , corresponding to in Definition 2.8, has a well-defined reaction network (for , , we restrict to positive integers), but an ill-defined kinetics taking the mass action form (we allow set to have elements that are negative). Thus, set corresponding to cannot be interpreted as a set of reaction rate constants, as opposed to set corresponding to (see also Example 2.2). Note also that , with .
2.3 Properties of kinetic functions
From Definition 2.3 it follows that a kinetic function has a structural property: cross-negative terms are absent. In this subsection, further properties of are defined: nonnegativity (absence of cross-negative effect), and a structural property called -factorability.
Let be given by , where , for , and If , , , then and are said to be nonnegative. Otherwise, and are said to be negative, and a cross-negative effect is said to exists for which such that and .
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 cross-negative 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 , one can prove that in two-dimensional 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).
Consider the following ODE system with polynomial RHS:
where , , and . Considering and , it follows that . Then:
System (2.5)–(2.6) induces a reaction system only in case (i). In particular, nonnegative ODE system (2.5)–(2.6) with 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 be given by , where , for , and Then term is said to be -factorable if , where is a polynomial function of . If , such that , then and ODE system are said to be -factorable. If it is true that , then and ODE system are said to be -factorable.
3 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 , 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 , while the second step is then to find a reaction system 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 often involve constraints defining simplicity of the system (e.g. see ), 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, , 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 , for a fixed kinetic function and a fixed set of complexes ( 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 , where 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.
Let be a kinetics. Consider the kinetic function given by , where and . Let us map to a reaction system with complexes and kinetics given by:
Reactant complexes, , are assumed to be uniquely obtainable from for , which is true in the case of the mass action kinetics.
Reaction is then constructed for each and , where new product complexes are given by , with being the sign function.
The new kinetics is then defined as , for , , where .
The induced reaction system is called the canonical reaction system, with 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.
3.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, , on the RHS (see Definition 2.7) . 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.
Let , , i.e. is a polynomial function. Consider the corresponing ODE system in the formal sum notation
where . Then, a transformation is called a kinetic transformation if the following conditions are satisfied:
, , , maps the polynomial function into a kinetic function for and .
Let be a fixed point of (3.2) that is mapped by to fixed point of the system of kinetic equations (2.1) with on its RHS. Let also the eigenvalues of the Jacobian matrix of , , denoted by , , be mapped to the eigenvalues of Jacobian of , , which are denoted by , . Then, for every such fixed point it must be true that , , and, if there are any additional eigenvalues , , they must satisfy .
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.
Let , . Let also be a continuous function, mapping set into , . Then, set is called a set of constraints.
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 -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 and higher 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 and lower ) 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.
Notation. Let , and , and suppose , and are each given a fixed ordering with indices being , , and , respectively, i.e. one can identify the ordered components of formal sums with components of Euclidean vectors. Let also the indices be denoted by , , 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 , , and .
3.2.1 Affine transformation
Consider applying an arbitrary nonsingular matrix on equation (3.3), resulting in:
where , and is a vector of new rate constants obtained from by rewriting the polynomial on the RHS of (3.4) into the mass action form. Then , mapping to , is called a centroaffine transformation. If is an orthogonal matrix, then is called an orthogonal transformation.
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 . 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 into new sets (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.
Let . If it is not possible that simultaneously is a kinetic function, and that a given set of constraints is satisfied, for all and for all , then it is said that 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 to be affinely kinetic, or nonkinetic, are generally difficult to obtain. Even in the simpler case , such conditions are complicated, and cannot be easily generalized for higher-dimensional systems and/or systems with higher degree of nonlinearily . In , 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.
If is nonkinetic under , given a set of constraints , for all orthogonal matrices and for all , then is also affinely nonkinetic, given , provided the following condition holds: , , for all diagonal and positive definite matrices , with .
By the singular value decomposition theorem, nonsingular matrices can be written as , where are orthogonal, and diagonal and positive definite. Cross-negative terms are invariant under transformation for all . If from Definition 3.3 is such that functions , , are invariant under all positive definite diagonal matrices , the statement of the theorem follows. ∎
3.2.2 X-factorable transformation
Consider multipling the RHS of equation (3.3) by a diagonal matrix , resulting in
Then , mapping to , is called an -factorable transformation. If is diagonal and its nonzero elements are
where , , then the transformation is denoted , and is said to be -factorable.
When is -factorable, i.e. , we write .
from Defnition 3.7 is a kinetic function, i.e. .
See . ∎
Functions and are not necessarily topologically equivalent due to two overlapping artefacts that can produce, so that is generally a nonkinetic transformation. Firstly, the fixed points of the former system can change the type and/or stability under , 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.
Let us consider the ODE system in two dimensions with RHS The following statements are true for all the fixed points of the two-dimensional system in under , , :
All the saddle fixed points are unconditionally invariant, i.e. saddle points of correspond to saddle points of .
A sufficient condition for stability of a fixed point to be invariant is:
A sufficient condition for the type of a fixed point to be invariant is:
Assume that the ODE system does not have fixed points on the axes of the phase space. Nevertheless, the two-dimensional system can have additional fixed points on the axes of the phase space, called boundary fixed points, denoted . The boundary fixed points can be either nodes or saddles, and the following statements are true:
If system is -factorable, then the origin is a fixed point, , with eigenvalues , , and the corresponding eigenvectors along the phase space axes.
For , , , , a boundary fixed point is a node if and only if
with the node being stable if , and unstable if , , . Otherwise, the fixed point is a saddle.
Without loss of generality, we consider two forms of the system (3.6) with :
where , so that system (3.7)–(3.8) is -factorable for , but only -factorable for . The results derived for an -factorable system hold when the system is -factorable, if the indices are swapped. By writing , , the Jacobian of (3.7)–(3.8), , is for given by
First, consider how fixed points of are affected by transformation . Denoting the Jacobian of two-dimensional system (3.3) by , and assuming the fixed points are not on the axes of the phase space (i.e. ), the Jacobians evaluated at are given by:
Comparing the trace, determinant and discriminant of and , we deduce (i)–(iii).
To prove (iv)–(v), we evaluate at the boundary fixed points of the form to get
If , then one of the boundary fixed points is , and the Jacobian becomes a diagonal matrix, so that condition (iv) holds. If , then in (3.9), and comparing the trace, determinant and discriminant of and , we deduce (v). ∎
Theorem 3.3 can be used to find conditions that an -factorable transformation given by is a kinetic transformation. While conditions (ii)–(iii) in Theorem 3.3 may be violated when is used, so that is a nonkinetic transformation, a composition of an affine transformation and an -factorable transformation, i.e. , may be kinetic. Furthermore, such a composite transformation may also be used to control the boundary fixed points introduced by . Finding an appropriate and 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 .
An alternative transformation, which is always kinetic, that also does not change the dimension of an ODE system is the time-parametrization transformation . However, while 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 , and introduces infinitely many boundary fixed points.
3.2.3 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 for any and
where , , and for all and . Assume further that the species set is partitioned, , , so that equations for species are kinetic, while those for species are nonkinetic. Consider the following total system, consisting of a degenerate system given by
which satisfies the initial condition (3.11) with for , and an adjointed system given by