Linear conjugacy of chemical kinetic systems

1.   Introduction

This paper presents a general computational solution to the problem of constructing linear conjugates of a chemical reaction network where each rate function is the product of a rate constant and an interaction function. We denote such a chemical kinetic system as a "rate constant-interaction function decomposable" (RID) kinetic system. Nearly all systems studied in Chemical Reaction Network Theory (CRNT) are RID kinetic systems, but recently "variable $k$" systems have been introduced in ^[1]. Furthermore, various kinetics such as weakly monotonic ones, are not explicitly required to have this form. Our approach is based on two new results:

1. The extension of the Johnston-Siegel Criterion for linear conjugacy (JSC) to the complex factorizable (CF) subset of RID kinetic systems, i.e., those whose interaction map $I_K: \Omega \rightarrow \mathbb{R}^\mathscr{R}$ factorizes via the space of complexes $\mathbb{R}^\mathscr{C}$: $I_K = I_k\circ \psi_K$ with $\psi_K: \Omega \rightarrow \mathbb{R}^\mathscr{C}$, with $\mathbb{R}_ > ^{\mathscr{S}}\subset \Omega \subset \mathbb{R}_\ge^\mathscr{S}$, as factor map and $I_k = diag(k)\circ \rho'$ with $\rho': \mathbb{R}^\mathscr{C} \rightarrow \mathbb{R}^\mathscr{R}$ assigning the value at a reactant complex to all its reactions (Theorem 4).

2. The dynamic equivalence of any non-complex factorizable (NF) RID kinetic system to a CF-system (Theorem 1).

An essential ingredient of the proofs of both results is the coincidence of the interaction maps of the kinetics considered. In the JSC extension (Theorem 4), the equality of the factor maps $\psi_K = \psi_K'$ (which is clearly equivalent to that of the interaction maps) is assumed. The CF-RM (Complex Factorizable by Reactant Multiples) transformation used to provide the dynamical equivalence in Theorem 1 is based on the concept of CF subsets of a reactant complex, which are defined as subsets of its reactions with the same interaction map. Determining the equality of functions (with infinite definition domains) may be computationally challenging, depending on their complexity and expression format. However, for a large subset of RID kinetic systems, which we call RID systems with interaction parameter maps (and denote with RIPK), the computational feasibility is ensured. Such systems are characterized by the existence of a map $P_K: \mathscr{R} \rightarrow \mathbb{R}^p$ such that $P_K = P_{K'}$ implies $I_K = I_{K'}$. The exponent $p$ is typically (but not always) a multiple of $m$ ( = number of species), and written as an $r \times p$ matrix. The interaction parameter map is easily seen as a generalization of the kinetic order matrix $F$ of power law kinetic systems.

Most RID kinetic systems, whose rate functions are specified explicitly, have interaction parameter maps, including all biochemical formalisms introduced to date. We discuss how the mixed integer linear programming (MILP) algorithms originally introduced for mass action kinetics (MAK) systems can be extended to RIP kinetic systems. We illustrate this and other results of the paper with an example of Hill-type kinetics (HTK), which was originally introduced as "Saturation Cooperativity Formalism" (SC Formalism) in ^[2].

The foundations for the study of dynamic equivalence in chemical kinetic systems were laid in the paper of Craciun and Pantea ^[3]. Important contributions to the theory in a more general context were previously provided by G. Farkas in ^[4]. The MILP-based computational approach to dynamic equivalence of MAK systems was pioneered by the group led by G. Szederk$\acute{e}$nyi and K. Hangos in Budapest, with further contributions from the lab of J. Banga in Vigo. Independently, M. Johnston and D. Siegel initiated the study of linear conjugacy, which led to the JSC for MAK systems. The three groups then collaborated in extending the MILP approach to linear conjugacy (a detailed discussion of the work up to 2013 can be found in ^[5]). Further developments included the extension to "Bio-CRNs" (whose rate functions are mass action functions divided by positive polynomials in the species variables) by G$\acute{a}$bor et al. ^[6] and to complex factorizable power law kinetic systems (denoted by PL-RDK) by Cortez et al. ^[7].

The paper is organized as follows: Section 2 collects the fundamentals of chemical reaction networks and kinetic systems required for the later sections. The central concept of "CF subsets of a reactant complex" and the method based on it are introduced in Section 3. The first main result (Theorem 1) is proved using the transformation. A Subspace Coincidence Theorem for the kinetic and stoichiometric subspaces (KSSC) of NF kinetic systems further illustrates the usefulness of CF-RM. Section 4 formulates the linear conjugacy problem for RID kinetic systems and extends the Johnston-Siegel Criterion (JSC) for linear conjugacy to complex factorizable RID systems. This is combined with the CF-RM method to provide the general computational solution to construct linear conjugates of any RID system. A running example (Examples 2–4), in Sections 3 and 4, further demonstrates the usefulness of the computational solution by deriving the existence of complex balanced equilibria of an NF power law kinetic system through construction of a weakly reversible, deficiency one PL-TIK system which is linear conjugate to the CF-RM transform. Section 5 focusses on the large subset of RID systems which have interaction parameter maps, for which the computational solution is always feasible. Details of the MILP-based algorithm are provided in Section 6. Section 7 illustrates the results of the paper using a reference system introduced in ^[2]. Conclusions and an outlook constitute Section 8. Tables of acronyms and frequently used symbols are provided in Supplementary Materials.

2.   Materials and method

We recall the necessary concepts of chemical reaction networks and the mathematical notation used throughout the paper adopted from the papers ^[7,8,9,10].

2.1.   Fundamentals of chemical reaction networks

We begin with the definition of a chemical reaction network.

Definition 1. A chemical reaction network is a triple $\mathscr{N} = (\mathscr{S}, \mathscr{C}, \mathscr{R})$ of three non-empty finite sets:

1. A set species $\mathscr{S}$,

2. A set $\mathscr{C}$ of complexes, which are non-negative integer linear combinations of the species, and

3. A set $\mathscr{R} \subseteq \mathscr{C} \times \mathscr{C}$ of reactions such that

● $(y, y) \notin \mathscr{R}$ for all $y \in \mathscr{C}$, and

● for each $y \in \mathscr{C}$, there exists a $y' \in \mathscr{C}$ such that $(y, y') \in \mathscr{R}$ or $(y', y) \in \mathscr{R}.$

We denote with $m$ the number of species, $n$ the number of complexes and $r$ the number of reactions in a CRN.

A complex is called monospecies if it consists of only one species, i.e., of the form $kX_i$, $k$ a non-negative integer and $X_i$ a species. It is called monomolecular if $k = 1$, and is identified with the zero complex for $k = 0$. A zero complex represents the "outside" of the system studied, from which chemicals can flow into the system at a constant rate and to which they can flow out at a linear rate (proportional to the concentration of the species). In biological systems, the "outside" also stands for the degradation of a species.

A chemical reaction network $(\mathscr{S}, \mathscr{C}, \mathscr{R})$ gives rise to a digraph with complexes as vertices and reactions as arcs. However, the digraph determines the triple uniquely only if an additional property is considered in the definition: $\mathscr{S} = \bigcup${ supp $y$ for $y \in \mathscr{C}\}$, i.e., each species appears in at least one complex. With this additional property, a CRN can be equivalently defined as follows.

Definition 2. A chemical reaction network is a digraph $(\mathscr{C}, \mathscr{R})$ where each vertex has positive degree and stoichiometry, i.e., there is a finite set $\mathscr{S}$ (whose elements are called species) such that $\mathscr{C}$ is a subset of $\mathbb{Z}^{\mathscr{S}}_{\geq}.$ Each vertex is called a complex and its coordinates in $\mathbb{Z}^{\mathscr{S}}_{\geq}$ are called stoichiometric coefficients. The arcs are called reactions.

Two useful maps are associated with each reaction:

Definition 3. The reactant map $\rho: \mathscr{R} \rightarrow \mathscr{C}$ maps a reaction to its reactant complex while the product map $\pi: \mathscr{R} \rightarrow \mathscr{C}$ maps it to its product complex. We denote $|\; \rho (\mathscr{\pi})\; |$ with $n_r$, i.e., the number of reactant complexes.

Connectivity concepts in Digraph Theory apply to CRNs, but have slightly differing names. A connected component is traditionally called a linkage class, denoted by $\mathscr{L}$, in CRNT. A subset of a linkage class where any two elements are connected by a directed path in each direction is known as a strong linkage class. If there is no reaction from a complex in the strong linkage class to a complex outside the same strong linkage class, then we have a terminal strong linkage class. We denote the number of linkage classes with $l$, that of the strong linkage classes with $sl$ and that of terminal strong linkage classes with $t$. Clearly, $sl \geq t \geq l.$

Many features of CRNs can be examined by working in terms of finite dimensional spaces $\mathbb{R}^{\mathscr{S}}, \mathbb{R}^{\mathscr{C}}, and \mathbb{R}^{\mathscr{R}},$ which are referred to as species space, complex space, and reaction space, respectively. We can view a complex $y \in \mathscr{C}$ as a vector in $\mathbb{R}^{\mathscr{C}}$ (called complex vector) by writing $y = \sum _{x \in \mathscr{S}} y_x x$, where $y_x$ is the stoichiometric coefficient of species $x$.

Definition 4. The reaction vectors of a CRN $(\mathscr{S}, \mathscr{C}, \mathscr{R})$ are the members of the set ${\{y'-y \in \mathbb{R}^{\mathscr{S}}\; |\; (y, y') \in \mathscr{R}\}}.$ The stoichiometric subspace $S$ of the CRN is the linear subspace of $\mathbb{R}^{\mathscr{S}}$ defined by

The rank of the CRN, $s$, is defined as $s = dim\; S.$

Definition 5. The incidence map $I_a: \mathbb{R}^{\mathscr{R}} \rightarrow \mathbb{R}^{\mathscr{C}}$ is defined as follows. For $f: \mathscr{R} \rightarrow \mathbb{R}$, then ${I_a(f)(v) = - f(a)}$ and $f(a)$ if $v = \rho(a)$ and $v = \pi(a)$, respectively, and are $0$ otherwise.

Equivalently, it maps the basis vector $\omega_a$ to $\omega_{v'} - \omega_v$ if $a: v \rightarrow v'$. It is clearly a linear map, and its matrix representation (with respect to the standard bases $\omega_a$, $\omega_{v}$) is called the incidence matrix, which can be described as

Let $I$ be the incidence matrix of the directed graph $D = (V, E)$. Then rank $I = n -l$, where $l$ is the number of connected components of $D$. A non-negative integer, called the deficiency, can be associated to each CRN. This number has been the center of many studies in CRNT due to its relevance in the dynamic behavior of the system. The deficiency of a CRN is the integer $\delta = n-l-s$. The reactant subspace $R$ is the linear space in $\mathbb{R}^\mathscr{S}$ generated by the reactant complexes. Its dimension, denoted by $q$, is called the reactant rank of the network. Meanwhile, the reactant deficiency $\delta_p$ is the difference between the number of reactant complexes and the reactant rank, i.e., $\delta_p = n_r -q$.

2.2.   Fundamentals of chemical kinetic systems

We now introduce the fundamentals of chemical kinetic systems. We begin with the general definitions of kinetics from ^[11]:

Definition 6. A kinetics for a CRN $\mathscr{N} = (\mathscr{S}, \mathscr{C}, \mathscr{R})$ is an assignment of a rate function $K_j: \Omega_K \rightarrow \mathbb{R}_\geq$ to each reaction $r_j \in \mathscr{R}$, where $\Omega_K$ is a set such that $\mathbb{R}^\mathscr{S}_ > \subseteq \Omega_K \subseteq \mathbb{R}^\mathscr{S}_\geq$, $c\wedge d \in \Omega_K$ whenever $c, d \in \Omega_K,$ and

A kinetics for a network $\mathscr{N}$ is denoted by ${K = (K_1, K_2, ..., K_r):\Omega_K \to {\mathbb{R}}^{\mathscr{R}}_{\geq}}$. A chemical kinetics is a kinetics $K$ satisfying the positivity condition: for each reaction $r_j:y\rightarrow y', K_j(c) > 0$ iff $\text{supp} y\subset\text{supp} c$. The pair $(\mathscr{N}, {K})$ is called the chemical kinetic system (CKS).

In the definition, $c \wedge d$ is the bivector of $c$ and $d$ in the exterior algebra of $\mathbb{R}^\mathscr{S}$. Once a kinetics is associated with a CRN, we can determine the rate at which the concentration of each species evolves at composition $c$.

Power-law kinetics is defined by an $r \times m$ matrix $F = [F_{ij}],$ called the kinetic order matrix, and vector $k \in \mathbb{R}^\mathscr{R}$, called the rate vector. In power-law formalism, the kinetic orders of the species concentrations are real numbers.

Definition 7. A kinetics $K: \mathbb{R}^\mathscr{S}_ > \rightarrow \mathbb{R}^\mathscr{R}$ is a power-law kinetics (PLK) if

with $k_i \in \mathbb{R}_ >$ and $F_{ij} \in \mathbb{R}.$

Definition 8. A chemical kinetics $K:\Omega \rightarrow \mathbb{R}^\mathscr{R}_\ge$ is complex factorizable (CF) if there is $k\in \mathbb{R}^\mathscr{R}_ >$ and a mapping $\psi_K: \Omega \rightarrow \mathbb{R}^{\mathscr{C}}$ such that $K = I_k\circ \psi_K$, where $I_k$ is the $\bf{k-}$interaction map defined by $I_k: \mathbb{R}^\mathscr{C} \rightarrow \mathbb{R}^\mathscr{R}$. The set of complex factorizable kinetics is denoted as $\mathscr{CFK(N)}$.

It can be deduced from the definition that if a chemical kinetics $K$ is complex factorizable, then its complex formation rate function $g = A_k\circ \psi_K$ and its species formation rate function (SFRF) $f = Y\circ A_k\circ \psi_K$. The $f(x) = \frac{dx}{dt}$ is the ODE or dynamical system of the {CKS}. A zero of $f$ is an element $c$ of $\mathbb{R}^\mathscr{S}$ such that $f(c) = 0$. A zero of $f$ is called an equilibrium (or steady state) of the ODE system. The SFRF contains three maps: map of complexes, Laplacian map, and factor map.

Definition 9. The map of complexes $Y: \mathbb{R}^\mathscr{C} \rightarrow \mathbb{R}^\mathscr{S}$ is defined by its values on the standard basis $\{\omega_y\}$, $y$ a non-zero complex: $Y(\omega_y) = y$ and extending it linearly to all elements of $\mathbb{R}^\mathscr{C}$. Its matrix, denoted with $Y$ (called the matrix of complexes), is an $m \times n$ matrix, its rows indexed by the species and its column by the complexes, with $y_{ij}$ being the stoichiometric coefficient of the $j^{th}$ complex in the $i^{th}$ species. In other words, the columns are the complexes written as column vectors.

Definition 10. The linear transformation $A_k: \mathbb{R}^\mathscr{C} \rightarrow \mathbb{R}^\mathscr{C}$ called Laplacian map is the mapping defined by $A_kx: = \sum_{(i, j)\in \mathscr{R}} k_{ij}x_i(\omega_j-\omega_i)$, where $x_i$ refers to the $i^{\text{th}}$ component of $x\in \mathbb{R}^\mathscr{C}$ relative to the standard basis. Its matrix representation is the $n\times n$ matrix such that

The label $k_{ji}$ is called the rate constant and is associated to the reaction $(j, i)\in \mathscr{R}$.

Definition 11. The factor map $\psi_K: \Omega \rightarrow \mathbb{R}^\mathscr{C}$ is defined as

Definition 12. A positive equilibrium or steady state $x$ is an element of $\mathbb{R}^\mathscr{S}_ >$ for which $f(x) = 0$. The set of positive equilibria of a chemical kinetic system is denoted by $E_+(\mathscr{N}, {K})$.

Two networks are said to be linearly conjugate if the solutions of their dynamic equations can be transformed into each other by a positive linear transformation ^[10,12].

Definition 13. Let $\Phi(x_0, t)$ and $\psi (x_0, t)$ be flows associated to kinetic systems $M$ and $M'$ respectively. $M$ and $M'$ are said to be linearly conjugate if there exists a bijective linear mapping $h : \mathbb{R}^n_{ > 0} \rightarrow \mathbb{R}^n_{ > 0}$ such that $h(\Phi(x_0, t)) = \psi(h(x_0, t))$ for all $x_0 \in \mathbb{R}^n_{ > 0}$.

Remark 1. In ^[13], it is shown that the bijection h in the previous definition corresponds to multiplication with a diagonal matrix with positive diagonal entries. The diagonal entries form the conjugacy vector $c$. More precisely, if $N$, $N'$ are the stoichiometric matrices and $K$, $K'$ are the kinetics of the systems $M$ and $M'$ respecively, then they are linearly conjugate if and only if $NK = \text{diag} (c)N'K'$.

Linear conjugacy is a generalization of the concept of dynamical equivalence.

Definition 14. Two kinetic systems are dynamically equivalent if the conjugacy vector $c = (1, \cdots, 1)$, i.e., if $NK = N'K'$.

In relation to linear conjugacy, if the mapping $h$ is trivial, $M$ and $M'$ are said to be dynamically equivalent ^[7].

2.3.   Rate constant-Interaction map Decomposable (RID) kinetics

To date, nearly all chemical kinetics studied in CRNT have constant rates, i.e. for each reaction $r$, the kinetic function $K_r: \Omega_K \rightarrow \mathbb{R}^\mathscr{R}$ can be written in the form $K_r(x) = k_r I_{K, r}(x)$, with a positive real number $k_r$ (called a rate constant) and an interaction map $I_{K, r}$. Recently however, G. Craciun and collaborators ^[1,14] have introduced variable $k$ systems, where the rates may vary between an upper and lower bound. Furthermore, there are kinetics sets such as the weakly monotonic kinetics studied in ^[15] or the span surjective kinetics introduced in ^[16] which do not explicitly require constant rates. The fractal kinetics studied primarily by physical chemists, e.g. Brouers ^[14] have rate values given by a function of exponential type. In view of this, we introduce the term Rate constant-Interaction map Decomposable (RID) kinetics for all chemical kinetics with constant rates and denote the set with RIDK.

In ^[8] (see also ^[16]), we introduce a special subset of $\mathscr{CFK(N)}$, which is the set of power law kinetics with reactant-determined kinetic orders, denoted by $\mathscr{PL-RDK (N)}$. A PLK system has a reactant-determined kinetic orders (of type PL-RDK) if for any two reactions $i, j$ with identical reactant complexes, the corresponding rows of kinetic orders in $V$ are identical, i.e., $v_{ik} = v_{jk}$ for $k = 1, 2, ..., m$.

We note also in ^[16] that $\mathscr{PL-RDK(N)}$ includes mass action kinetics (MAK) and coincides with the set of GMAK systems recently introduced by ^[17] if the vertices map $y: \mathscr{C} \rightarrow R^m$ of the GMAK system is injective. They also constitute the subset of power law systems for which various authors claimed that their results "hold for the complexes with real coefficients" are valid.

Another important property of a complex factorizable kinetics is "factor span surjectivity":

Definition 15. Let $f:V\rightarrow W$ be a map between finite dimensional vector spaces $V$ and $W$. $f$ is span surjective if and only if $\text{span}(\text{Im} f) = W$.

In ^[16], it is shown that $f$ is span surjective if and only if its coordination functions are linearly independent.

Definition 16. A complex factorizable kinetics $K$ is factor span surjective if its factor map $\psi_K$ is span surjective. $\mathscr{FSK(N)}$ denotes the set of factor span surjective kinetics on a network $\mathscr{N}$.

We characterized in ^[16] a factor span surjective PL-RDK system.

Proposition 1. A PL-RDK system is factor span surjective if and only if no rows corresponding in the kinetics order matrix $F$ corresponding to different reactant complexes coincide (i.e. $\rho(r)\ne \rho(r')\Rightarrow F_r\ne F_{r'}$).

We recall the definition of the $m \times n$ matrix $\tilde{Y}$ from ^[17]: for a reactant complex, the column of $\tilde{Y}$ is the transpose of the kinetic order matrix row of the complex' reaction, otherwise (i.e. for a terminal point), the column is 0.

The $\bf{T}-$matrix of a PL-RDK system is formed by truncating away the columns of the terminal points in $\tilde{Y}$, obtaining an $m \times n_r$ matrix. The corresponding linear map $T : \mathbb{R}^{\rho(\mathscr{R})} \rightarrow \mathbb{R}^{\mathscr{R}}$ maps $\omega_\rho(r)$ to $(F_r)^T$. The subspace $\tilde{R} : = \text{ Im }T = \big \langle (F_r)^T \big \rangle$ is called the kinetic reactant subspace and $\tilde{q} = \dim \tilde{R}$ is called the kinetic reactant rank of the system.

Let $e^1, e^2, ..., e^\ell \in \left\{ 0, 1\right\}^{n}$ be the characteristic vectors of the sets $\mathscr{C}^{1}$, $\mathscr{C}^{2}$, ..., $\mathscr{C}^{\ell}$, respectively, where $\mathscr{C}^{i}$ is the set of complexes in linkage class $\mathscr{L}^{i}$. That is, for all $j \in \mathscr{C}$ and $i = 1, \dots, \ell$, we have $e^i_j = 1$ if $j \in \mathscr{C}^i$, and 0 otherwise. Let $L = \left[e^1, e^2, ..., e^\ell \right]$. Define the $\bf{\hat{T}}-$matrix, an $(m+\ell)\times n_{r}$ block matrix, by

where $L_{pr}$ is the truncated matrix $L$ (i.e., non-reactant rows are left out).

If the non-inflow columns (i.e., columns of the complexes associated to non-inflow reactions) of ${T}-$matrix corresponding to each linkage class are linearly independent and if its column rank is maximal, then the chemical kinetics is said to be $\bf{\hat{T}}-$rank maximal (to type PL-TIK).

3.   CF Transformation of NF Kinetics

The CF-RM (Complex Factorization by Reactant Multiples) method developed from a proposal by C. Pantea in December 2017 for such a transformation for power law kinetics. The key idea is, at an NF branching point, i.e. a complex which is the reactant of reactions (called its branching reactions) with non-proportional interaction maps, to transform reactions by introducing new reactants while conserving the reaction vectors, thus leaving the stoichiometric subspace invariant. CF-RM refines the approach by ensuring that the reactant subspace also remains invariant and that a minimum number of reactions is transformed. The essential underlying concept of CF-RM is that of CF-subsets of the set of reactions of a reactant complex. The concept is also the basis for the construction of CF-decompositions of a RID kinetic system.

3.1.   CF-subsets of the reaction set of a reactant complex

For a reactant complex $y$ of a network $\mathscr{N}$, $\mathscr{R}(y)$ denotes its set of (branching) reactions, i.e., $\rho^{-1}(y)$ where $\rho:\mathscr{R} \rightarrow \mathscr{C}$ is the reactant map. The $n_r$ reaction sets $\mathscr{R}(y)$ of reactant complexes partition the set of reactions $\mathscr{R}$ and hence induce a decomposition of $\mathscr{N}$.

Definition 17. Two reactions $r$, $r' \in \mathscr{R}(y)$ are CF-equivalent for K if their interaction functions coincide, i.e., $I_{K, r} = I_{K, r'}$ or, equivalently, if their kinetic functions $K_r$ and $K_r'$ are proportional (by a positive constant). The equivalence classes are the CF-subsets (for $K$) of the reactant complex $y$.

Definition 18. If $N_R(y)$ is the number of CF-subsets of $y$, then ${1 \le N_R(y) \le \mid \rho^{-1}(y)\mid}$. The reactant complex is a CF-node if $N_R(y) = 1$, and an NF-node otherwise. It is a maximally NF-node if ${N_R(y) = \mid\rho^{-1}(y)\mid > 1}$.

Definition 19. The number N$_\boldsymbol{R}$ of CF subsets of a CRN is the sum of $N_R(y)$ over all reactant complexes.

Clearly, $N_R \ge n_r$ and the kinetics $K$ is CF if and only if $N_R = n_r$, or equivalently all reactant complexes are CF-nodes for $K$.

Example 1. For a power law kinetic system, the CF-subsets of a reactant complex are the subsets of branching reactions with identical rows in the kinetic order matrix. To show this, we recall that the interaction map of a PLK system is $= x^F$ and hence the claim is $x^{l(r)} = x^{l(r')} \Rightarrow l(r) = l(r')$. The "$\leq$" is evident, for the converse, let $e_i$ be the positive vector with e (the exponential number) as its ith coordinate and 1's otherwise. Since $\log x^{l(r)} = l(r)\log x$, the value of the $\log$ at $e_i =$ the ith kinetic order, which proves the claim.

Example 2. (Running Example - Part 1) In ^[18], a power law kinetic system for R. Schmitz's pre-industrial carbon cycle model was introduced. The system (depicted in Figure 1) with 6 complexes (representing carbon pools) and 13 reactions (indicating mass transfer) is weakly reversible and has zero deficiency.

The system's kinetic order matrix is given by:

The column to the right of F lists the rate constants of the corresponding reactions. The kinetic order matrix reveals that the system has 3 NF nodes (reactant complexes): $M_1$, $M_2$ and $M_3$. The following table lists their CF-subsets:

Hence, $\mathscr{R}$ is partitioned into 9 CF-subsets, i.e., $N_R = 9$.

We also note that since the CF-subsets of a reactant complex partition its reaction set, and the reaction sets of reactant complexes partition the set of reactions, that the CF-subsets determine a decomposition.

We recall from ^[19] that a subset $\mathscr{R}'$ of $\mathscr{R}$ defines a subnetwork $\mathscr{N}' = (\mathscr{S}', \mathscr{C}', \mathscr{R}')$ with $\mathscr{C}'$ consisting of the complexes occurring in reaction of $\mathscr{R}'$ and $\mathscr{S}'$ consisting of the species occurring in complexes in $\mathscr{C}'$. A CRN decomposition $\mathscr{N} = \mathscr{N}_1 \cup... \cup \mathscr{N}_k$ consists of the subnetworks $\{ \mathscr{N}_i \}$ induced by a partition $\{ \mathscr{R}_i \}$ of $\mathscr{R}$. We use the model presented in ^[20] to illustrate the concepts introduced above.

Definition 20. The CF-subsets of a RID kinetic system partition the reaction set and induce the CFS decomposition of the system.

The CFS decomposition consists of $N_R$ subnetworks, whereby $n_r \le N_R \le r$.

3.2.   CFM decompositions of a RID kinetic system

In this Section, we introduce useful coarsenings of the CFS-decomposition of a RID kinetic system.

For each NF node $y$, we choose an ordering of its CF-subsets $\mathscr{R}_1(y)$, $\mathscr{R}_2 (y)$, ..., $\mathscr{R}_{NR}(y)$ according to decreasing number of reactions in the CF-subset. We define $\mathscr{R}_i : = \cup_{y\in \rho(R)} \mathscr{R}_i(y)$ where $i = 1, ..., \max_{y\in \rho(R)} N_R(y)$ and $\mathscr{R}_i(y') = \phi$ if $N_R(y') < i$.

We can now introduce the concept of a maximal CF-subsystem (CFM) of a RID kinetic system:

Definition 21. A maximal CF-subsystem $(\mathscr{N}_{mcf}, {K})$ of a RID kinetic system $(\mathscr{N}, {K})$ is induced the union of the reaction sets of all CF-nodes and a CF-subset with the maximal number of reactions from each NF-node, i.e., the union $\mathscr{R}_{mcf}$ of $\{ \rho^{-1}(y)\mid y\ is \ CF-node\}$ and $\mathscr{R}_1$.

Clearly, there may be several maximal CF-subsystems in a RID kinetic system, but the number of reactions in each of them is the same, and we denote this with $r_{mcf}$. Note that since $\mid \mathscr{R}_i\mid \ge \mid \mathscr{R}_j\mid$ if $i < j$, then $r_{mcf}\ge \mid \mathscr{R}_i\mid$ for all $i$.

Definition 22. A CFM decomposition is induced by the reaction set partition $\{\mathscr{R}_{mcf}, \mathscr{R}_2, ..., \mathscr{R}_k\}$, with $k = \max_{y\in \rho(R)} N_R(y)$.

A CFM-decomposition is clearly a coarsening of the CFS-decomposition. It is the decomposition into CF-subsystems with the least number of subnetworks.

3.3.   CF-RM Transformation of an NF kinetic system: the generic case

We first introduce the concept of a CF-transformation of an NF kinetic system:

Definition 23. A CF kinetic system $(\mathscr{N}^*, K^*)$ is a CF-transform of an NF system $(\mathscr{N}, {K})$, where $\mathscr{N} = (\mathscr{S}, \mathscr{C}, \mathscr{R})$, $\mathscr{N}^{*} = (\mathscr{S}^{*}, \mathscr{C}^{*}, \mathscr{R}^{*})$ and $N$, $N^{*}$ as their respective stoichiometric matrices, if and only if $\mathscr{S}^{*} = \mathscr{S}$, $N^{*} = N$, and $K^{*} = K$.

$N^{*}K^{*} = NK$ implies that a CF-transform is dynamically equivalent to the original NF system. Moreover, the stoichiometric subspaces coincide, i.e., $S^* = S$.

Our first main result is the following Theorem:

Theorem 1. Any NF system $(\mathscr{N}, {K})$ is dynamically equivalent to a CF system $(\mathscr{N}^{*}, {K}^{*})$ via a CF-transformation.

Proof. We construct the CF-transformation nodewise, i.e., we specify how to transform each NF-node $y$ into $N_R(y)$ CF-nodes. Let $\mathscr{R}_1(y), ..., \mathscr{R}_k(y)$ (where $k = N_R(y)$) be the CF subsets of $y$. We leave $\mathscr{R}_1(y)$ unchanged. We choose a complex $y_2$ such that $y + y_2$ is not contained in $\rho(\mathscr{R})$. All reactions in $\mathscr{R}_2(y)$ are transformed "catalytically", i.e., $r_i: y \rightarrow z_i$ is replaced by $r_i'´: y + y_2 \rightarrow z_i + y_2$. The reaction vector is unchanged. For the reactions in $\mathscr{R}_3(y)$, choose a complex $y_3$ such that $y + y_3$ is not in $\rho(\mathscr{R}) \cup \{ y + y_2\}$ and proceed as in $\mathscr{R}_2(y)$. After $N_R(y) -1$ steps, we have completed the transformation for $y$. After the transformation of all NF nodes, we have a CF-transform as claimed.

There is clearly a multitude of ways to carry out CF-transformations, and a good principle is to minimize the changes needed as well as keep further network components invariant under the necessary changes. In this spirit, the specific goals of the CF-RM method are:

● minimize the number of reactions to be changed and

● leave the reactant subspace invariant, i.e., $R^{*} = R$.

The first goal is achieved by choosing, for each NF node, a CF subset with the maximal number of reactions, as the subset to be left unchanged. The second goal is accomplished by selecting the "catalytic" complexes used as multiples of the reactant complex (as expressed in the acronym CF-RM).

The CF-RM method proceeds as follows:

● Determine the reactant set $\rho(\mathscr{R})$ (see Algorithm 1 lines 1–4).

● A CF-node is left unchanged (see Algorithm 1 lines 5–21).

● At an NF-node, select a CF-subset with the maximal number of reactions. Note that there may be several. This CF-subset is left unchanged (this step minimizes the number of t-reactions overall and may see Algorithm 1 lines 22–29).

● For each of the remaining $N_R(y) - 1$ CF-subets, choose successively a multiple of $y$ which is not among the current set of reactants, i.e., those of the original networks left unchanged and the already selected new reactants. Various procedures are possible for this selection of a new reactant; the essential condition is that it is different from those in the current reactant set. After each choice, the current set must be updated. For each Non-reactant Determined Kinetics (NDK) reactant complex $y$, $N_R(y) - 1$ new reactants are constructed (see Algorithm 1 lines 30–37).

● Since the last expression is also true for a CF-node, the total number of new reactants $= \sum (N_R(y) - 1)$ with the sum taken over all reactants. This number $= \sum N_R(y) - \sum 1 = \sum N_R(y) - n_r = N_R - n_r$. Under CF-RM the number of CF-subsets NR of the original system is also the number of reactants of the transformed system, since the latter is equal to $n_r + N_R - n_r = N_R$.

Remark 2. If an NF system has at least one NF-node with more than 1 CF-subset with the maximal number of reactions, then several transforms can be generated, which might have some differing network properties. It is possible to define an additional procedure for which CF-subset to choose and leave unchanged.

Remark 3. As mentioned above, various procedures can be defined to select a new reactant. One possible procedure is the following:

● Determine the set of multiples of $y$ among the current reactants.

● If the set is empty, set $m_y = 1$.

● If the set is non-empty, determine the maximum multiple $y' = max_y y$. Set $m_y = max_y$.

● The new reactant is $y + m_y y$.

Instead of repeating the reactant set check for every CF-subset of $y$, one could further optimize by ordering the CF-subsets to be changed, doing the above for the first, and then use $y + (m_y + i - 1)y$ for the $i = 2, ..., N_R-1$.

Table 1 presents the key network numbers of a CF-RM transform in equations or inequalities involving only network numbers of the original NF network. Thus, the relationships are of predictive character.

Remark 4. The addition of complexes to both sides of a reaction is similar to the technique used by M. Johnston for translating mass action systems to generalized mass action systems in ^[21].

In the next proposition, we provide a proof of a Table 1 entry which is not straightforward.

Proposition 2. $i)$ $l^{*} = l^{*}_r + l^{*}_b + l$, where $l^{*}_r =$ number of new linkage classes generated by new reactants and $l^{*}_b =$ number of new linkage classes due to link-breaking.

$ii)$ $l^{*} - l^{*}_b \le (N_R - n_r) + l$.

Proof. For $i)$, the equation expresses the partitioning into 3 subsets. For $ii)$, a new reactant adds at most 1 linkage class (none if it coincides with an old product complex or at least one of the new product complexes in its linkage class coincides with an old complex).

The "link-breaking" effect of CF-RM is shown in the CRN in Figure 2: if $R_{(i-1)}: X_1 \rightarrow X_i$ for $i = 2, ..., 5$, $R_5: X_4 \rightarrow X_6$, $R_6: X_5 \rightarrow X_7$ and $X_1$ NF with CF-subsets $\{R_1, R_2\}$ and $\{R_5, R_6\}$, then $\delta^{*} = 10 - 4 - 6 = 0 = \delta$.

One notes however that three key network numbers of $\mathscr{N}^{*}$ have question marks: the number of complexes $n^{*}$, the deficiency $\delta^{*}$ and the number of terminal strong linkage classes $t^*$. Indeed, for many networks, the deficiency increases under CF-RM, but, as the following Proposition shows, for certain network classes, it decreases.

Proposition 3. Let $d$ be an integer $\ge 2$. Let $\mathscr{N}_d$ be the CRN with species $X_1$, $X_2$ and the following reactions:

$R_1: X_1 \rightarrow 2X_1$

$R_i: X_1 \rightarrow 2iX_1 + X_2$ for $i = 2, ..., d$

$R_{d+i -1}: X_1 \rightarrow (2i - 1)X_1 + X_2 \rightarrow X_1 + (2i-1)X_2: R_{2d+i -2}$ for $i = 2, ..., d$

Let $X_1$ be an NF node with CF-subsets $\{R_1, \cdots, R_d\}$ and $\{R_{d+1}, \cdots, R_{2d-1}\}$.

Then, $\delta - \delta^{*} = d - 1$.

Proof. The new reactions are: $2X_1 \rightarrow 2iX_1 + X_2$ for $i = 2, ..., d$. The remaining reactant complexes are all non-branching, thus RDK and unchanged. Hence there is no new complex, while there are $d -1$ new linkage classes due to the "link-breaking" effect, i.e., $n^{*} = n, l^{*} = 1 + (d -1) = d \Rightarrow \delta^{*} = n - d - 2$. Therefore, $\delta - \delta^{*} = (n - 1 - 2) - (n - 2 - d) = d - 1$.

In the next section, we present a special variant of CF-RM where these network numbers can be better estimated.

3.4.   CF-RM$_\boldsymbol{+}$: a "choosier" CF-RM variant

CF-RM$_+$ is a variant of CF-RM which uses additional criteria in the selection of the new reactant multiples. All other steps are identical with the generic CF-RM method, i.e., a CF-RM$_+$ transform is also a CF- transform.

CF-RM$_+$ chooses the reactant multiple so that

a) the new reactant differs from all existing complexes, and

b) all the new product complexes in the CF-subset also differ from all existing complexes.

There are of course various ways of ensuring that conditions a) and b) are fulfilled and we leave it to the first consequence of transforming via CF-RM$_+$, which is a more predictable change in deficiency.

Proposition 4. For a CF-RM$_+$ transform $\mathscr{N}^{*}$, $\delta^{*}\ge \delta$.

Proof. For any CRN, $n = n_r + t_p$, where $t_p$ is the number of terminal points. In an CF-RM$_+$ transform, in each subset to be changed, there is one new reactant complex and exactly $x$ new terminal points. The number of reactions to be changed in the CF-subset is also pertained by $x$. Since all terminal point of the original network are conserved (with no coincidence), we obtain $n^{*} - n = (N_R - n_r) + (r - r_{mcf})$. On the other hand, $l^{*} - l = l^{*}_r + l^{*}_b$. For any CF-RM$_+$ transform, $l^{*}_b \le r - r_{mcf}$ (a link-break is created per new reactant–whether it leads to a new linkage class or not depends on specific network properties). This implies that $l^{*} - l \le (N_R - n_r) +(r - r_{mcf})$. Hence, $\delta^{*}- \delta = (n^{*}- l^{*}) - (n - l) = (n^{*} - n) - (l^{*} - l) \ge (N_R - n_r) + (r - r_{mcf})- (N_R - n_r) - (r - r_{mcf}) = 0$.

Remark 5. The monomolecular system from Figure 2 shows this lower bound is sharp.

Besides the change in deficiency, the change in the number of terminal strong linkage classes is difficult to predict under the generic CF-RM transformation. Recall that $t$ has two components, i.e., $t = t_p + t_c$, which are the number of terminal points and the number of cycle terminal classes. Under CF-RM$_+$, the relationships for its components can be predicted and together provide an expression for the change in $t$ as shown in the following Proposition:

Proposition 5. For a CF-RM$_+$ transform $\mathscr{N}^{*}$, we have:

i) $t_p^{*} - t_p = r - r_{mcf}$

ii) $t_c^{*} - t_c \le 0$

iii) $t^{*} - t \le r - r_{mcf}$

Proof. $i)$ was already shown (and used) in the previous Section. For $ii)$ note that a reversible pair of reactions can be broken up into two irreversible reactions under CF-RM$_+$. On the other hand, no new cycles can emerge since there is no coincidence of new complexes with existing ones. $iii)$ follows by adding $i)$ and $ii)$.

Corollary 1. For a CF-RM$_+$ transform $\mathscr{N}^{*}, n^{*} = n + (N_R - n_r) + (r - r_{mcf})$.

Proof. In the identity $n^{*} - n = (n_r^{*} - n_r) + (t_p^{*} - t_p)$, we substitute $N_R$ for $n_r^{*}$ and use Proposition 5.i.

Example 3. (Running Example - Part 2) To apply CF-RM to Schmitz's carbon cycle model, we replace $R_3$, $R_4$, and $R_7$ with the following reactions:

Each of the new reactions forms a linkage class of $\mathscr{N}^{*}$, with the remaining original 10 reactions of $\mathscr{N}$ forming the fourth one depicted in Figure 3:

Table 2 presents the network numbers of the $\mathscr{N}^{*}$.

The network is $t$-minimal, but clearly not weakly reversible (in fact, it is point-terminal). Note that it is also a CF-RM$_+$ transform.

The $T$ matrix of the CF system $(\mathscr{N}^{*}, K^*)$ is given by:

3.5.   A Subspace Coincidence Theorem for NF kinetic systems

In this Section, we present an initial application of CF-RM transformation by deriving a Subspace Coincidence Theorem for NF systems.

Arceo et al. ^[16] generalized the Subspace Coincidence Theorem of Feinberg and Horn ^[9] from MAK systems to CF systems as follows:

Theorem 2. For a complex factorizable system on a network $\mathscr{N}.$

1) If $t-l > \delta$, then $K\ne S$.

1') If $0 < t-l = \delta$, and a positive steady state exists, then $K\ne S$. In fact $dim \ S-dim \ K\ge t-l-\delta+1$ if the system is also factor span surjective.

2) If $t-l = 0$ (i.e., $\mathscr{N}$ is $t-$minimal), then $K = S$.

3) If $0 < t-l < \delta$ or $= \delta$ and a positive steady state does not exist, then it is rate constant dependent whether $K = S$ or not.

We first note that for any CF-RM transform $\mathscr{N}^{*}$, we not only have coincidence of stoichiometric subspaces $S = S^{*}$ but also coincidence of the kinetic subspaces $K = K^{*}$ (due to the dynamic equivalence, $f = f^*$, implying $\text{Im} f = \text{Im} f^{*}$ and $\text{span} (\text{Im} f) = \text{span} (\text{Im} f^{*})$.

Our approach is to identify properties for an NF system so that its CF-RM$_+$ transform satisfies the conditions of the Theorem above. Our first step is to extend the kinetics concept of factor span surjectivity, which is currently defined only for CF systems, to any RID kinetic system.

A CF-subset $\mathscr{R}_i$ is characterized by the common interaction map $I_K(\mathscr{R}_ i)$ of the kinetics of its reactions. This leads to the following definition:

Definition 24. A RID kinetics is interaction span surjective if and only if the set $\{I_K(\mathscr{R}_ i)\}$ of its CF-subset interaction maps is linearly independent.

The following Proposition shows that "interaction span surjectivity" is the correct extension of the factor span surjectivity concept.

Proposition 6. If $(\mathscr{N}, K)$ is interaction span surjective, then its CF-transform $(\mathscr{N}^{*}, K^{*})$ is also factor span surjective.

Proof. Since $(\mathscr{N}^{*}, K^{*})$ is CF, $N_R^{*} = n_r^{*}$. On the other hand, the latter is equal to $N_R$. Hence the set of interaction maps of $\mathscr{N}$ and $\mathscr{N}^{*}$ coincide. For a CF system, since $I_K = \rho' \cdot \psi_K$, it is clear that linear independence of both sets are equivalent.

As a second step, we identify the network properties of the NF-system $(\mathscr{N}, K)$ such that the properties needed to apply the various statements of the Theorem to $(\mathscr{N}^{*}, {K}^{*})$ are ensured.

We first state two Lemmas.

Lemma 1. If $\mathscr{N}$ is SRD, then $\mathscr{N}^{*}$ is also SRD.

Proof. $n^{*}_r = N_R \ge n_r \ge s = s^{*}$.

The second Lemma is a general relationship between TBD and SRD networks derived from a (submitted) manuscript by Farinas et al. entitled "Species subsets and embedded networks of S-systems":

Lemma 2. Let $\mathscr{N}$ be a chemical reaction network.

i) A network with deficiency-bounded terminality (TBD) has sufficient reaction diversity.

ii) If the network is point terminal, then the converse also holds, i.e., $TBD\le SRD$ (or equivalently $TND\ge LRD$).

We can now state and prove a Subspace Coincidence Theorem for NF-systems:

Theorem 3. Let $(\mathscr{N}, {K})$ be an NF RIDK system.

1) If $N_R < s$, then $K \ne S$.

If the system is also intersection span surjective, then either

2) $\mathscr{N}$ is $t-$minimal and $r - r_{mcf} = N_R - n_r$, implies $K = S$; or

3) $\mathscr{N}$ is TBD and point terminal, implies that $K = S$ is rate-constant dependent.

Proof. $1)$ $n^{*}_r = N_R < s = s'$ means that $\mathscr{N}^{*}$ is LRD. By Lemma 2 $(i)$, it follows that it is also TND, and by (1) of the KSSC in ^[12], $K = K^{*} \ne S^{*} = S$.

$2)$ In order to apply (2) of the Arceo et al. KSSC ^[12], we need to show that there is a CF-RM transform such that $\mathscr{N}^{*}$ is $t-$minimal, or $t^{*} - l^{*} = 0$. We calculate this difference for an CF-RM$_+$ transform as follows: $t^{*} - l^{*} = t^{*}_p + t^{*}_c - l^{*} = r - r_{mcf} + t_p + t^{*}_c- l^{* }$(by Proposition 5) $= r - r_{mcf} + t_p + t^{*}_c - ((N_R-n_r) + t_p + t_c)$ (based on the properties of CF). After canceling terms, we obtain $0 \le t^{*} - l^{*} = t^{*}_c - t_c \le 0$ (by Proposition 5), implying the claim.

$3)$ $\mathscr{N}$ is point terminal $\Rightarrow \mathscr{N}^{*}$ is point terminal (by Proposition 5.ⅱ), hence after Lemma 1 and Lemma 2 $(ii)$, $\mathscr{N}^{*}$ is also TBD, and (3) of Arceo et al. Theorem can be applied, implying $K = K^{*} = S^{*} = S$ is rate constant dependent.

Remark 6. Since both the stoichiometric and reactant subspaces of an NF system and its CF-RM transform coincide, the underlying networks have the same R and S class introduced in ^[12]. This implies that a Theorem for the coincidence of kinetic and reactant subspaces of NF systems analogous to that for CF-systems derived in ^[12] can also be stated and proved.

4.   Linear conjugacy of RID kinetic systems

In this Section, we present a solution to the problem of finding linear conjugates of any RID kinetic system. After extending the Johnston-Siegel Criterion (JSC) for linear conjugacy to CF systems, we can generate linear conjugates for any RID kinetic system by applying the JSC to any CF-RM transform of the given system. We also discuss some computational challenges regarding the solution approach.

4.1.   The Johnston-Siegel Criterion for linear conjugacy (JSC) of CF kinetic systems

Theorem 4. Consider two CF systems $\left(\mathscr{N}, {K}\right)$ and $\left(\mathscr{N'}, {K'}\right)$ with $\mathscr{N} = \left(\mathscr{S}, \mathscr{C}, \mathscr{R} \right)$ and $\mathscr{N}' = \left(\mathscr{S'}, \mathscr{C'}, \mathscr{R'} \right)$. Let $Y = Y'$ be the matrix of complexes for both networks. Suppose further that the factor maps coincide, i.e., $\psi_K = \psi_{K}'$. Let $A_b$ be a Laplacian with the same structure as that of $\left(\mathscr{N'}, \mathscr{K'}\right)$ and $c$, a positive vector in $\mathbb{R}^m$ such that $Y\cdot A_k = C\cdot Y\cdot A_b$, where $C = \text{diag} (c)$. Then $\mathscr{N}$ is linearly conjugate to $\mathscr{N'}$ with the Laplacian $A_k' = A_b \cdot \text{diag} (\psi_{K'(c)})$.

Proof. Let $\varphi \left({{x_o}, t} \right)$ be the solution of the system of ODE $\dot{x} = f\left(x \right) = Y\cdot{A_k}\cdot{\psi _K}$ associated to the reaction network $\mathscr N$.

Consider the linear map $h\left(x \right) = {C^{ - 1}}x$ where $C = diag(c)$, $c \in \mathbb{R}_{ > 0}^n$.

Let $\tilde\varphi \left({{y_0}, t} \right) = {C^{ - 1}}\varphi \left({{x_0}, t} \right)$ so that $\varphi \left({{x_0}, t} \right) = C\tilde\varphi \left({{y_0}, t} \right)$. It follows that

Now,

where $D = diag (e)$ and ${e_j} = \left\{ {\begin{array}{*{20}{c}} {{c^{{F._j}}}{\text{, if complex }}j{\text{ is a reactant of some reaction }}k{\text{ }}} \\ {1{\text{, otherwise}}} \end{array}} \right.$

So, $\tilde \varphi ' \left({{y_0}, t} \right) = Y\cdot{A_b}\cdot D\cdot{\psi _K}\left({\tilde \varphi \left({{y_0}, t} \right)} \right)$. Clearly, $\tilde \varphi \left({{y_0}, t} \right)$ is a solution of the system $\dot{x} = Y\cdot{A_b}\cdot D\cdot{\psi _K}$ corresponding to the reaction network $\mathscr N'$. We have that $h\left({\varphi \left({{x_0}, t} \right)} \right) = \tilde \varphi \left({h\left({{x_0}} \right), t} \right)$ for all ${x_0} \in \mathbb{R}_{ > 0}^n$ and $t\ge 0$ where $y_0 = h(x_0)$ since ${y_0} = \tilde \varphi \left({{y_0}, t} \right) = {C^{ - 1}}{y_0} = \varphi \left({{y_0}, t} \right)$. It follows that networks $\mathscr N$ and $\mathscr N'$ are linearly conjugate.

4.2.   A solution to the linear conjugacy problem of RID kinetic systems

A solution approach to the linear conjugacy problem of RID kinetic systems is clearly to first transform the system if necessary (i.e., if it is an NF system) via CF-RM to a CF system and then apply the Johnston-Siegel Criterion to generate linearly conjugate systems. The second step could be done using MILP algorithms based on the JSC, once these are extended to appropriate CF systems (cf. Section 6).

Example 4. (Running Example - Part 3) In ^[18], Fortun et al. derived a Deficiency Zero Theorem for a class of NF power law kinetic systems and applied it to a subsystem of the Schmitz's carbon cycle model to establish the existence of positive equilibria for the subsystem. The authors then used a "Lifting Theorem" of ^[19] to show the existence of corresponding positive equilibria for the whole system. Here, we provide an alternative approach for this result by using the MILP algorithm of ^[7], a special case of the MILP algorithm introduced in Section 6, to construct a weakly reversible PL-TIK system, which is linearly conjugate to the CF-transform of Schmitz's model discussed previously. The results of ^[22] show that this weakly reversible system has positive equilibria, and hence so does its linear conjugate, the Schmitz's carbon cycle model.

The sparse linear conjugate of $(\mathscr{N}^*, K^*)$ was obtained using the MILP algorithm, described in ^[7]. The algorithm seeks to generate linearly conjugate realizations for a class of power-law kinetic systems, i.e., PL-RDK. Prior to the implementation of the algorithm, the map of complexes $Y$, the Laplacian map $A_k$, and kinetic order matrix $F$ are required to be set first. The matrix $F$ was given in the preceding section. The following are the associated matrices $Y$ and $A_k$ of the system.

where $C_1: M_1$, $C_2: M_2$, $C_3: M_3$, $C_4: 2M_1$, $C_5: M_5+M_1$, $C_6: 2M_2$, $C_7: M_1+M_2$, $C_8: M_4$, $C_9: 2M_3$, $C_{10}: M_1+M_3$, $C_{11}: M_5$, and $C_{12}: M_6$.

Additionally, the parameters were set as follows: $\epsilon = 0.001$ and ${u_{ij} = 20, i, j = 1, 2, ..., 12, i\ne j}$. Using MATLAB R2018b, the linearly conjugate weakly reversible sparse realization ($\tilde{\mathscr{N}}, \tilde{ K}$) was obtained with the corresponding Laplacian map $A_k^{sparse}$.

The linear conjugacy constants are $c_1 = 2.28$, $c_2 = 1.14$, $c_3 = 1.14$, $c_4 = 1.14$, $c_5 = 4.56$, and $c_6 = 4.56$. Furthermore, the associated system of ODEs is given below:

The network $\tilde{\mathscr{N}}$ and its numbers are shown in the following Figure 4 and Table 3:

The $\hat{T}$ matrix of the system is given by:

One readily computes that it has maximal rank, 9, and hence $(\tilde{\mathscr{N}}, \tilde{K})$ is a PL-TIK system. Since each of the linkage classes has zero deficiency, according to the Deficiency Zero Theorem for PL-TIK systems (Theorem 5 and Corollary 6, ^[22]), each subsystem possesses positive equilibria. It then follows from Theorem 4 of ^[22] that the whole system also has positive equilibria. Hence, the linearly conjugate system $(\mathscr{N}, K)$ also has positive equilibria, which are necessarily complex balanced since the system has zero deficiency. The graphs of the individual trajectories of ($\mathscr{N}^ *, K^*$) and $(\tilde{\mathscr{N}}, \tilde{K})$ are depicted in Figure 5.

There are however several challenges with this "solution in principle": It may be difficult to compute the CF subsets of a RID kinetic system, which form the basis of the CF-RM method, as it involves determining if interaction functions (for an infinite number of domain values) are equal. This clarity depends on how explicit and complex the functional expressions are. Similarly, applying the JSC to a CF system, one needs to establish the equality of the factor maps, which is equivalent to the difficulty with interaction functions cited above.

In the next section, we identify a large subset of RID kinetics, where the solution approach can be applied in general.

5.   Linear conjugacy of RIP kinetic systems

This section introduces the large subset of RID kinetics with interaction parameter maps (RIPK). The subset includes power law kinetics (PLK), Hill-type kinetics (HTK)–originally called "Saturation-Cooperativity" (SC) Formalism ^[2], and other published biochemical kinetics such as linlog ^[23] and loglin kinetics ^[24]. We extend the $T$ matrix concept of ^[22] to complex factorizable RIP kinetics (denoted by RIP-CFK) and obtain a computationally feasible form of the JSC for this kinetics set, which leads to executable solutions of the linear conjugacy problem.

5.1.   RIP kinetics: RID kinetics with interaction parameter map

Definition 25. A set $\mathscr{K} \in$ RIDK is said to be of type "RID kinetics with interaction parameter maps" if there is a family of maps $\left\{ p_{\mathscr{K}}: \mathscr{R} \rightarrow R^{m1} \times... \times R^{mk} \mid K \in \mathscr{K}\right\}$ such that

i) $p_K(r) = p_K(r') \Rightarrow I_K(x)_r = I_K(x)_{r'}$ for all $x$ in $\Omega$ and

ii)$< /p > < p > ● p_K = p_{K'} \Rightarrow I_K(x) = I_{K'}(x)$ for all $x$ in $\Omega$

Example 5. PLK with the family of kinetic order matrices, i.e., $p_\mathscr{K}(r) = F_r$, (kinetic order row vector or interaction), is the primary example. Since $I_K(x) = x^F$, the properties i) and ii) are straightforward.

Example 6. Hill-type kinetics (HTK)–originally called "Saturation-Cooperativity (SC) Formalism" in 2007 by Sorribas et al. ^[2]. We recall the definition of ^[16]:

Definition 26. Hill-type kinetics (HTK) is defined as follows:

with $c\in \mathbb{R}^n_\ge$ (defined by continuity at the boundary), $k_j\in \mathbb{R}_ >$, $d_j\in \mathbb{R}^n_\ge$ and $v_j\in \mathbb{R}^n$ for $j = 1, ..., m$. Note that the $v_j$ have to be nonegative.

The family of interaction parameter maps is given by $P_K: \mathscr{R} \rightarrow \mathbb{R}^m \times \mathbb{R}^m$ with $p_K(r) = (v_1, ..., v_m, d_1, ..., d_m)$, where we leave out the index $j$.

5.2.   CF-RM for RIP-NFK and the JSC for RIP-CFK

Since under CF-RM, there is a bijection $\eta : \mathscr{R} \rightarrow \mathscr{R}^*$, if $(\mathscr{N}, {K})$ is an NF RIP kinetic system, then $(\mathscr{N}^*, {K})$ is a CF RIP kinetic system with the interaction parameter map $p_{K}^* (\eta (r)) : = p_K(r)$.

We denote the set of all complex factorizable kinetics with interaction parameter maps with RIP-CFK.

For an interaction parameter map $p_K: \mathscr{R} \rightarrow R_{m1} \times... \times R_{mk}$, we write $p = m1 +...+ mk$. It is now easy to formally introduce the $T$ matrix of a RIP-CFK kinetics:

Definition 27. The T matrix of a RIP-CFK kinetics $K$ is the $p \times n_r$ matrix whose jth column is $p_K(r)^T$, where $\rho(r) = j$. The $\hat{T}$ matrix is the $(p + l) x$ given by adjoining the characteristic functions of the linkage classes as rows to the T matrix. The rank of the $\hat{T}$ matrix is denoted by $\hat{q}$.

We have the following useful Proposition:

Proposition 7. Let $(\mathscr{N}, {K})$ and $(\mathscr{N}', K')$ be RIP-CFK systems. If $T = T'$, then $\psi_K = \psi_{K'}$.

Proof. $T = T' \Rightarrow p_K = p_{K'}$ for all $K, K'$ of the same type $\Rightarrow I_K = I_{K'}$ (by definition of interaction parameter map) $\Leftrightarrow \psi_K = \psi_{K'}$ (since the maps differ only with the reactions map). Hence, RIP-CF kinetics, it suffices to check a finite set of vectors to establish the coincidence of the factor maps. This allows the extension of the JSC-based MILP algorithms for PL-RDK systems to RIP-CFK systems. Since the CF-RM transform of a RIP-NFK system is clearly a RIP-CFK system, we obtain a general computational solution for the linear conjugacy of RIP kinetic systems.

Remark 7. The set $\{\mathscr{K} \in RIPK \mid \rho(r) = \rho(r') \Rightarrow p_K(r) = p_K(r')\}$ may, in general, be a proper subset of RIP-CFK. This may result in computing a smaller set of linear conjugates as when the whole set RIP-CFK is used. This is a small price one pays for ensuring the computational feasibility. There are, however, various RIP kinetics for which the converse $\psi_K(x) = \psi_{K'}(x) \Rightarrow p_K(r) = p_K(r')$ also holds, so that the corresponding sets are equal. Examples are PLK and $PYK_h$ (the set of poly-PL kinetics with h summands), which form a covering of PYK (cf. a manuscript in preparation by Talabis et al. entitled "A Weak Reversibility Theorem for poly-PL kinetics and the replicator equation").

In ^[16], we introduced the notations PL-RDK and HT-RDKD for the subsets of PLK and HTK respectively, which satisfy $\rho(r) = \rho(r') \Rightarrow p_K(r) = p_K(r')$. For any other subset A of RIPK, we will denote $\{\mathscr{K} \in RIPK \mid \rho(r) = \rho(r') \Rightarrow p_K(r) = p_K(r')\}$ with A-RDP (kinetics with reactant-determined parameter maps). This notation is consistent with earlier ones since the corresponding letters there indicate the specific parameter maps, too.

6.   Extension of MILP algorithms to RIP-CFK systems

Cortez et al. ^[7] extended the MILP algorithm developed by Johnston et al. ^[25] to find linearly conjugate networks of PL-RDK systems. Aside from linear conjugacy, other desirable properties can be incorporated in the algorithm such as weak reversibility and minimal deficiency (e.g. deficiency zero). In this study, we focus on extending the algorithm to find linearly conjugates of RIP-CFK kinetic systems.

6.1.   Key components of the MILP algorithm

The algorithm considers two CF systems: the original system $(\mathscr{N}, {K})$ with $\mathscr{N} = (\mathscr{S}, \mathscr{C}, \mathscr{R})$ and the target system $(\mathscr{N}', {K}')$ with $\mathscr{N}' = (\mathscr{S}', \mathscr{C}', \mathscr{R}')$. The algorithm determines the corresponding network structure of the target system satisfying the linear conjugacy property. The two networks $\mathscr{N}$ and $\mathscr{N}'$ have the same set of species and complexes. As a consequence, their corresponding molecularity matrices and the coefficient maps coincide. The algorithm requires that $\mathscr{R}$ and ${K}$ be known while $\mathscr{R}'$ and ${K}'$ are to be obtained. The following are needed to be ascertained prior to the MILP implementation:

● molecularity matrix $Y\in \mathbb{R}^{m\times n}_{\ge 0}$;

● matrix $M = Y \cdot A_k$, where $A_k$ is the Laplacian map;

● parameter $\epsilon > 0$, that is set to be sufficiently small; and

● parameter $u_{ij} > 0$, where $i, j = 1, ..., m$, $i \ne j$.

Remark 8. Note that $\epsilon$ and $u$ are introduced to ensure the correct structure of the linearly conjugate realization.

6.2.   MILP algorithm to CF systems

The MILP algorithm finds a sparse linearly conjugate realization of the original network $\mathscr{N}$. A sparse realization contains the minimum number of reactions, hence the associated objective function of the MILP model is

There are two sets of constraints in the model which indicate the linear conjugacy condition and desired structure of the network.

Table 4 shows the description of the variables used in the model. Constraint (6.2) imposes the linear conjugacy specification while constraint (6.3) ensures that the target network $\mathscr{N}'$ has the correct structure. A dense linearly conjugate network can also be determined by considering the maximization problem analog.

The optimal solution (if it exists) of the MILP would yield the matrix $A_b$ with the same structure as $\mathscr{N}'$, and the conjugacy constant vector $c$. The Laplacian map $A_k'$ of the target network is computed as:

where $D = diag(e)$ and ${e_j} = \left\{ {\begin{array}{*{20}{c}} {{c^{F \cdot j}}, {\text{if complex $j$ is a reactant of some reaction }}k}\\ {1, \text{otherwise}} \end{array}} \right.$

6.3.   MILP algorithm to RIP-CFK systems

It is important to note that the algorithm developed by Cortez et al. ^[7] is only applicable to CF systems (e.g. PL-RDK). For NF systems, the MILP cannot be immediately utilized to generate linearly conjugate realizations. It is necessary to transform it into a CF system through the CF-RM algorithm described in Section 5. This framework is applicable to RIPK systems which include both the power law kinetics (PLK) and Hill-type kinetics (HTK). The computation of the matrix $A_b$ and linear conjugacy vector $c$ is the same for both systems. The process of finding linearly conjugate realizations differs only in the derivation of corresponding sytem of ODEs wherein the respective kinetic order matrix/interaction parameter matrix is incorporated accordingly. Additionally, to obtain a proper form of the rational terms in the target HTK system, a linear scaling of the variable of rational term must be carried out, that is the variable must be multiplied by its corresponding linear conjugacy constant. This approach is similar to the approach of ^[26] to linear conjugacy of bio-CRNs.

7.   Application to Hill-type kinetic system

In ^[16] and ^[8], we introduced CRN representations of GMA systems–as defined in Biochemical Systems Theory (BST)–by means of the biochemical maps usually used to define them. These representations are actually independent of the power law kinetics assigned to the reactions from BST and we will use them for other RID kinetic systems too, as illustrated in the following examples.

In the following, after a brief review of Hill-type kinetics, we consider a reference metabolic system of ^[2]. We apply the MILP algorithm to Hill-type kinetics and compare the set of linear conjugates with those of power law kinetics on the same chemical reaction network.

7.1.   Review of Hill-type kinetics

The set of Hill-type kinetics was introduced in 2007 by Sorribas et al. ^[2] under the name of "Saturation-Cooperativity Formalism" (SC-Formalism). This framework generalizes the well-known Michaelis-Menten and Hill functions in one variable. The term "Hill-type kinetics" (HTK) was introduced in 2013 in the paper of Wiuf and Feliu ^[10]. In ^[16], it was shown that a Hill-type kinetics can be written as follows:

Given

● $e_j:\Omega_K\rightarrow \mathbb{R}^\mathscr{S}_ >$ with $e_j(x) = (x^{\rho(j)_1}_1, ..., x^{\rho(j)_m}_m), r_j\in\mathscr{R}$

● $d_j: \mathbb{R}^\mathscr{S}_ > \rightarrow \mathbb{R}^\mathscr{S}_ >$ with $d_j(x) = x+D_j$, $r_j\in \mathscr{R}$

● $m: \mathbb{R}^\mathscr{S}_ > \rightarrow \mathbb{R}_ >$ with $m(x) = \Pi x_i, i = 1, ..., m.$

then $I_H = I_1/I_2$, with $I_1$ a PLK interaction map with kinetic order matrix $F$ and $(I_2 (x))_ j = m \cdot d_j \cdot e_j (x)$. Furthermore, the dissociation vectors $d_j$ (s. Definition 22) were organized in an $r \times m$ matrix called the "dissociation matrix" and the set of complex factorizable Hill-type kinetics was denoted by HT-RDKD (Hill-type with reactant-determined kinetic and dissociation), expressing the fact that it is the pre-image of the interaction parameter map given by the kinetic order and dissociation matrices.

Remark 9. The method for determining linear conjugates for Bio-CRNs in ^[6] is applicable to HTK if the exponents are non-negative integers.

7.2.   The reference system with Hill-type kinetics

Now, we apply the integrated algorithm to a particular biological system. Specifically, we consider a metabolic network with one positive feedforward and a negative feedback (see Figure 6) taken from the published work of ^[2].

The corresponding embedded representation of the metabolic network, with $X_5$ as an independent variable, is as follows:

We apply the MILP algorithm on the SC Formalism approximation by ^[2] of the reference model depicted in Figure 6. Using the framework, the corresponding system of ODEs for the reference model is given as:

where $V_1 = 8$, $V_2 = 84.2175$, $V_3 = 8$, $V_4 = 115.341$, $V_5 = 8$, and $V_6 = 8$. The interaction parameter matrix (containing the kinetic orders and dissociation constants) for the given system is:

with $n_{21} = 1$, $n_{23} = -0.8429$, $n_{32} = 1$, $n_{41} = 2.9460$, $n_{42} = 3$, $n_{53} = 1$, $n_{64} = 1$, $k_{21} = 0.6705$, $k_{41} = 0.8581$, $k_{42} = 44.7121$, $k_{53} = 1$, and $k_{64} = 1$.

Using the parameter values $u_{ij} = 20, i, j = 1, 2, ..., 9$ for $i\ne j$ and $\epsilon = 0.1$ and considering the same matrices $Y$ and $M$, the sparse linearly conjugate network of the Hill-type system is

with the corresponding system of ODEs

where $\bar{V}_1 = 0.8$, $\bar{V}_2 = 12.0921$, $\bar{V}_3 = 8$, $\bar{V}_4 = 10185531.88$, $\bar{V}_5 = 8$, and $\bar{V}_6 = 8$.

The linearly conjugate dense realization was also obtained. The structure of the network is given as:

The conjugacy constants of the derived network are: $c_1 = 2.9555$, $c_2 = 9.9140$, $c_3 = 0.4$, and $c_4 = 10$. Using these constants and the computed $A_b$, we obtained the corresponding Kirchhoff matrix for the network:

The associated ODEs for the dense realization is

with $\bar{W}_1 = 2.7068$, $\bar{W}_2 = 54.3516$, $\bar{W}_3 = 127.9649$, $\bar{W}_4 = 7.0009$, $\bar{W}_5 = 52.067$, $\bar{W}_6 = 0.9914$, $\bar{W}_7 = 8$, $\bar{W}_8 = 2372.74$, $\bar{W}_9 = 273674.2$, and $\bar{W}_{10} = 8$. The kinetic orders and dissociation constants are $n_{21} = n_{31} = 1$, $n_{23} = n_{33} = -0.8429$, $n_{42} = n_{52} = n_{62} = 1$, $n_{73} = 1$, $n_{81} = n_{91} = 2.9460$, $n_{82} = n_{92} = 3$, and $n_{104} = 1$, $k_{21} = k_{31} = 0.6705$, $k_{23} = k_{33} = 3.9065$, $k_{42} = k_{52} = k_{62} = 1$, $k_{73} = 1$, $k_{81} = k_{91} = 0.8581$, $k_{82} = k_{92} = 44.7121$, and $k_{104} = 1$.

The linearly conjugate sparse network has also 6 reactions which is equal to the number of reactions of the derived linearly conjugate sparse system with power-law kinetics. Whereas, the dense realization of the SC model has 10 reactions. The graphs of the individual trajectories of the original Hill-type system and the linearly conjugate systems are depicted in S1a-S1d and Figures S1e-S1h, respectively.

8.   Conclusion

Different networks could generate the same set of ODEs making them dynamically equivalent. In the past few years, various authors have pioneered the use of MILP algorithms for determining linear conjugacy between MAK systems ^[5,13,25,27], between rational functions systems ^[26], between GMAK systems ^[28] and between PL-RDK systems ^[7]. In the work of ^[7], they extended the JSC for linear conjugacy from MAK systems to PL-RDK systems. It is limited to power law kinetic systems with branching reactant complexes that have identical kinetic orders. In this study, we further extended the algorithm for branching reactant complexes with different kinetic orders.

We summarize below main results presented in this paper:

1. We showed that any non-complex factorizable (NF) RID kinetic system can be dynamically equivalent to a CF system via CF-transformation (Theorem 1).

2. We further illustrated the usefulness of CF-RMA through the extended proof of Subspace Coincidence Theorem for the kinetic and stoichiometric subspaces (KSSC) of NF kinetic systems.

3. We extended the JSC for linear conjugacy to the CF subset of RID kinetic systems, i.e., those whose interaction map $I_K: \Omega \rightarrow \mathbb{R}^\mathscr{R}$ factorizes via the space of complexes $\mathbb{R}^\mathscr{C}$: $I_K = I_k\circ \psi_K$ with $\psi_K: \Omega \rightarrow \mathbb{R}^\mathscr{C}$ as factor map and $I_k = diag(k)\circ \rho'$ with $\rho': \mathbb{R}^\mathscr{C} \rightarrow \mathbb{R}^\mathscr{R}$ assigning the value at a reactant complex to all its reactions (Theorem 4).

4. We demonstrated (with running examples: Examples 2 - 4) that linear conjugacy can be generated for any RID kinetic systems by applying the JSC to any NF kinetic system that are transformed to CF kinetic system. The extended JSC for linear conjugacy to CF-RID systems is combined with the CF-RM method to provide the general computational solution to construct linear conjugates of any RID system.

5. For a large subset of RID kinetic systems RIPK, which have interaction parameter maps, we illustrated how the proposed approach of this paper can also be applied and that the computational solution is always feasible. We presented an example of HTK which was also known as SC Formalism.

Acknowledgments

We thank Casian Pantea for presenting the idea of transforming any power law kinetic system to a dynamically equivalent reactant-determined system, which was the basis for the development of the CF-RM method. ARL held research fellowships from De La Salle University and would like to acknowledge the support of De La Salle University's Research Coordination Office.

Conflict of interest

We have no conflicts of interest to disclose.

Supplementary