Linear conjugacy of chemical kinetic systems

Two networks are said to be linearly conjugate if the solution of their dynamic equations can be transformed into each other by a positive linear transformation. The study on dynamical equivalence in chemical kinetic systems was initiated by Craciun and Pantea in 2008 and eventually led to the Johnston-Siegel Criterion for linear conjugacy (JSC). Several studies have applied Mixed Integer Linear Programming (MILP) approach to generate linear conjugates of MAK (mass action kinetic) systems, Bio-CRNs (which is a subset of hill-type kinetic systems when the network is restricted to digraphs), and PL-RDK (complex factorizable power law kinetic) systems. In this study, we present a general computational solution to construct linear conjugates of any"rate constant-interaction function decomposable"(RID) chemical kinetic systems, wherein each of its rate function is the product of a rate constant and an interaction function. We generate an extension of the JSC to the complex factorizable (CF) subset of RID kinetic systems and show that any non-complex factorizable (NF) RID kinetic system can be dynamically equivalent to a CF system via transformation. We show 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.


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 IK : Ω → R R factorizes via the space of complexes R C : IK = I k •ψK with ψK : Ω → R C , with R S > ⊂ Ω ⊂ R S ≥ , as factor map and I k = diag(k)•ρ with ρ : R C → R 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 CFsystem (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 ψK = ψ 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 PK : R → R p such that PK = P K implies IK = I K .The exponent p is typically (but not always) a multiple of m (= number of species), and written as an r × 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é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á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.

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 -10 ).

Fundamentals of chemical reaction networks
We begin with the definition of a chemical reaction network.Definition 1.A chemical reaction network is a triple N = (S , C , R) of three non-empty finite sets: (1) A set species S , (2) A set C of complexes, which are non-negative integer linear combinations of the species, and • (y, y) / ∈ R for all y ∈ C , and • for each y ∈ C , there exists a y ∈ C such that (y, y ) ∈ R or (y , y) ∈ 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 kXi, k a non-negative integer and Xi 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 (S , C , 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: S = { supp y for y ∈ 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 (C , R) where each vertex has positive degree and stoichiometry, i.e., there is a finite set S (whose elements are called species) such that C is a subset of Z S ≥ .Each vertex is called a complex and its coordinates in Z S ≥ are called stoichiometric coefficients.The arcs are called reactions.
Two useful maps are associated with each reaction: Definition 3. The reactant map ρ : R → C maps a reaction to its reactant complex while the product map π : R → C maps it to its product complex.We denote | ρ(π) | with nr, 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 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 ≥ t ≥ l.
Many features of CRNs can be examined by working in terms of finite dimensional spaces R S , R C , andR R , which are referred to as species space, complex space, and reaction space, respectively.We can view a complex y ∈ C as a vector in R C (called complex vector ) by writing y = x∈S yxx, where yx is the stoichiometric coefficient of species x.The rank of the CRN, s, is defined as s = dim S. Definition 5.The incidence map Ia : R R → R C is defined as follows.For f : R → R, then , respectively, and are 0 otherwise.
Equivalently, it maps the basis vector ωa to ω v − ωv if a : v → v .It is clearly a linear map, and its matrix representation (with respect to the standard bases ωa, ω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 δ = n − l − s.The reactant subspace R is the linear space in R S generated by the reactant complexes.Its dimension, denoted by q, is called the reactant rank of the network.Meanwhile, the reactant deficiency δp is the difference between the number of reactant complexes and the reactant rank, i.e., δp = nr − q.

Fundamentals of chemical kinetic systems
We now introduce the fundamentals of chemical kinetic systems.We begin with the general definitions of kinetics from (11 ): A chemical kinetics is a kinetics K satisfying the positivity condition: for each reaction rj : y → y , Kj(c) > 0 iff supp y ⊂ supp c.The pair (N , K) is called the chemical kinetic system (CKS).
In the definition, c ∧ d is the bivector of c and d in the exterior algebra of R 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 × m matrix F = [Fij], called the kinetic order matrix, and vector k ∈ R R , called the rate vector.In power-law formalism, the kinetic orders of the species concentrations are real numbers.
with ki ∈ R> and Fij ∈ R.
and a mapping ψK : Ω → R C such that K = I k • ψK , where I k is the k−interaction map defined by I k : R C → R R .The set of complex factorizable kinetics is denoted as C F K (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 • ψK and its species formation rate function dt is the ODE or dynamical system of the CKS.A zero of f is an element c of R 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 : R C → R S is defined by its values on the standard basis {ωy} , y a non-zero complex: Y (ωy) = y and extending it linearly to all elements of R C .Its matrix, denoted with Y (called the matrix of complexes), is an m × n matrix, its rows indexed by the species and its column by the complexes, with yij 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 : R C → R C called Laplacian map is the mapping defined by A k x := (i,j)∈R kijxi(ωj − ωi), where xi refers to the i th component of x ∈ R C relative to the standard basis.Its matrix representation is the n × n matrix such that The label kji is called the rate constant and is associated to the reaction (j, i) ∈ R.
Definition 11.The factor map ψK : Ω → R C is defined as Definition 12.A positive equilibrium or steady state x is an element of R S > for which f (x) = 0.The set of positive equilibria of a chemical kinetic system is denoted by E+(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 Φ(x0, t) and ψ(x0, 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 : 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 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, In relation to linear conjugacy, if the mapping h is trivial, M and M are said to be dynamically equivalent (7 ) .

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 Kr : ΩK → R R can be written in the form Kr(x) = krIK,r(x), with a positive real number kr (called a rate constant) and an interaction map IK,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 C F K (N ), which is the set of power law kinetics with reactant-determined kinetic orders, denoted by 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 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 : C → 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 → W be a map between finite dimensional vector spaces V and W . f is span surjective if and only if span(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 ψK is span surjective.F S K (N ) denotes the set of factor span surjective kinetics on a network 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.ρ(r) = ρ(r ) ⇒ Fr = F r ).
We recall the definition of the m × n matrix Ỹ from (17 ): for a reactant complex, the column of Ỹ 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 T−matrix of a PL-RDK system is formed by truncating away the columns of the terminal points in Ỹ , obtaining an m × nr matrix.The corresponding linear map T : R ρ(R) → R R maps ωρ(r) to (Fr) T .The subspace R := Im T = (Fr) T is called the kinetic reactant subspace and q = dim R is called the kinetic reactant rank of the system.Let e 1 , e 2 , ..., e ∈ {0, 1} n be the characteristic vectors of the sets C 1 ,C 2 ,...,C , respectively, where C i is the set of complexes in linkage class L i .That is, for all j ∈ C and i = 1, . . ., , we have e i j = 1 if j ∈ C i , and 0 otherwise.Let L = e 1 , e 2 , ..., e .Define the T−matrix, an (m + ) × nr block matrix, by where Lpr 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 T−rank maximal (to type PL-TIK).

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.

CF-subsets of the reaction set of a reactant complex
For a reactant complex y of a network N , R(y) denotes its set of (branching) reactions, i.e., ρ −1 (y) where ρ : R → C is the reactant map.The nr reaction sets R(y) of reactant complexes partition the set of reactions R and hence induce a decomposition of N .Definition 17.Two reactions r, r ∈ R(y) are CF-equivalent for K if their interaction functions coincide, i.e., IK,r = I K,r or, equivalently, if their kinetic functions Kr and K r are proportional (by a positive constant).The equivalence classes are the CF-subsets (for K) of the reactant complex y.Clearly, NR ≥ nr and the kinetics K is CF if and only if NR = nr, 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 ) ⇒ l(r) = l(r ) .The "≤" is evident, for the converse, let ei 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 ei = 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 preindustrial 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): M1, M2 and M3.The following table lists their CF-subsets: Hence, R is partitioned into 9 CF-subsets, i.e., NR = 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 R of R defines a subnetwork N = (S , C , R ) with C consisting of the complexes occurring in reaction of R and S consisting of the species occurring in complexes in C .A CRN decomposition N = N1 ∪ ... ∪ N k consists of the subnetworks {Ni} induced by a partition {Ri} of 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 NR subnetworks, whereby nr ≤ NR ≤ r.

CFM decompositions of a RID kinetic system
In this Section, we introduce useful coarsenings of the CFS-decomposition of a RID kinetic system.
We can now introduce the concept of a maximal CF-subsystem (CFM) of a RID kinetic system: Definition 21.A maximal CF-subsystem (N mcf , K) of a RID kinetic system (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 R mcf of {ρ −1 (y) | y is CF − node} and R1.
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 | Ri |≥| Rj | if i < j, then r mcf ≥| Ri | for all i.Definition 22.A CFM decomposition is induced by the reaction set partition {R mcf , R2, ..., R k }, with k = max y∈ρ(R) NR(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.

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 (N * , K * ) is a CF-transform of an NF system (N , K), where N = (S , C , R), N * = (S * , C * , R * ) and N , N * as their respective stoichiometric matrices, if and only if S * = S , N * = N , and K * = K.N * K * = N K 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 (N , K) is dynamically equivalent to a CF system (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 NR(y) CF-nodes.Let R1(y), ..., R k (y) (where k = NR(y)) be the CF subsets of y.We leave R1(y) unchanged.We choose a complex y2 such that y + y2 is not contained in ρ(R).All reactions in R2(y) are transformed "catalytically", i.e., ri : y → zi is replaced by r i : y + y2 → zi + y2.The reaction vector is unchanged.For the reactions in R3(y), choose a complex y3 such that y + y3 is not in ρ(R) ∪ {y + y2} and proceed as in R2(y).After NR(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).
• 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 NR(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, NR(y) − 1 new reactants are constructed (see Algorithm 30-37).• Since the last expression is also true for a CF-node, the total number of new reactants = (NR(y) − 1) with the sum taken over all reactants.This number = NR(y) − 1 = NR(y)−nr = NR −nr.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 nr+NR−nr = NR.
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 my = 1.
• If the set is non-empty, determine the maximum multiple y = maxyy.Set my = maxy.
• The new reactant is y + myy.
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 + (my + i − 1)y for the i = 2, ..., NR − 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.INPUT2: column sum of ρ (from OUTPUT2) 7: OUTPUT3: identify the branching complexes return complex y is a branching reactant complex return complex y is a non-branching reactant complex INPUT3: kinetic order of the identified branching complex (from OUTPUT3) 17: OUTPUT4: determine whether the branching complex is RDK or NDK

18:
if all kinetic order associated to the identified branching complex are all equal then Let N R be the number of distinct kinetic order representation for each z.note: the reactions of this RDK-subset is left unchanged.

33:
check the reactant a in ρ(R)

35:
Let m a be the coefficient of a in ρ(R).

36:
OUTPUT9: transform the reactions in R b as such that the new reactant is a + m a a = (m a + 1)a 37: OUTPUT10: Update ρ(R) (from OUTPUT9) 38: REPEAT Procedure CF Transformation (for the remaining distinct kinetic order values) 39: REPEAT Procedure Generate RDK subsets for every NDK complexes (for the remaining NDK node) 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  One notes however that three key network numbers of N * have question marks: the number of complexes n * , the deficiency δ * 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 ≥ 2. Let N d be the CRN with species X1, X2 and the following reactions: Proof.The new reactions are: 2X1 → 2iX1 + X2 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 In the next section, we present a special variant of CF-RM where these network numbers can be better estimated.

CF-RM + : 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.Proof.For any CRN, n = nr+tp, where tp 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 = (NR − nr) + (r − r mcf ).On the other hand, l * − l = l * r + l * b .For any CF-RM+ transform, l * b ≤ 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 ≤ (NR − nr) + (r − r mcf ).Hence, 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 = tp + tc, 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 N * , we have: 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).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 (N * , K * ) is given by:

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 N .
3) If 0 < t − l < δ or = δ 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 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 Im f = Im f * and span(Im f ) = span(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 Ri is characterized by the common interaction map IK (Ri) 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 {IK (Ri)} 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 (N , K) is interaction span surjective, then its CF-transform (N * , K * ) is also factor span surjective.
On the other hand, the latter is equal to NR. Hence the set of interaction maps of N and N * coincide.For a CF system, since IK = ρ • ψ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 (N , K) such that the properties needed to apply the various statements of the Theorem to (N * , K * ) are ensured.
We first state two Lemmas.
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 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., T BD ≤ SRD (or equivalently T N D ≥ LRD).
We can now state and prove a Subspace Coincidence Theorem for NF-systems: Theorem 3. Let (N , K) be an NF RIDK system.
If the system is also intersection span surjective, then either 2) N is t−minimal and r − r mcf = NR − nr, implies K = S; or 3) N is TBD and point terminal, implies that K = S is rate-constant dependent.

Proof.
1) n * r = NR < s = s means that N * is LRD.By Lemma 2 (i), it follows that it is also TND, and by (1) of the KSSC in ( 12), K = K * = 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 N * is t−minimal, or t * − l * = 0. We calculate this difference for an CF-RM+ transform as follows: 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.

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.Let Y = Y be the matrix of complexes for both networks.Suppose further that the factor maps coincide, i.e., ψK = ψ K .Let A b be a Laplacian with the same structure as that of (N , K ) and c, a positive vector in R m such that Y Proof.Let ϕ (xo, t) be the solution of the system of ODE ẋ where D = diag(e) and ej = c F. j , if complex j is a reactant of some reaction k and t ≥ 0 where y0 = h(x0) since y0 = φ (y0, t) = C −1 y0 = ϕ (y0, t).It follows that networks N and N are linearly conjugate.

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).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 (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.Additionally, the parameters were set as follows: = 0.001 and uij = 20, i, j = 1, 2, ..., 12, i = j.Using MATLAB R2018b, the linearly conjugate weakly reversible sparse realization ( Ñ , K) was obtained with the corresponding Laplacian map The linear conjugacy constants are c1 = 2.28, c2 = 1.14, c3 = 1.14, c4 = 1.14, c5 = 4.56, and c6 = 4.56.Furthermore, the associated system of ODEs is given below: The network Ñ and its numbers are shown in the following Figure 4 and Table 3: The T matrix of the system is given by: One readily computes that it has maximal rank, 9, and hence ( Ñ , 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 (N , K) also has positive equilibria, which are necessarily complex balanced since the system has zero deficiency.The graphs of the individual trajectories of (N * , K * ) and ( Ñ , 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.

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 K ∈ RIDK is said to be of type "RID kinetics with interaction parameter maps" if there is a family of maps p K : R → R m1 × ... × R mk | K ∈ K such that i) pK (r) = pK (r ) ⇒ IK (x)r = IK (x) r for all x in Ω and ii) pK = p K ⇒ IK (x) = I K (x) for all x in Ω Example 5. PLK with the family of kinetic order matrices, i.e., p K (r) = Fr, (kinetic order row vector or interaction), is the primary example.Since IK (x) = x F , the properties i) and ii) are straightforward.
The family of interaction parameter maps is given by PK : R → R m × R m with pK (r) = (v1, ..., vm, d1, ..., dm), where we leave out the index j.

CF-RM for RIP-NFK and the JSC for RIP-CFK
Since under CF-RM, there is a bijection η : R → R * , if (N , K) is an NF RIP kinetic system, then (N * , K) is a CF RIP kinetic system with the interaction parameter map p * K (η(r)) := pK (r).
We denote the set of all complex factorizable kinetics with interaction parameter maps with RIP-CFK.
For an interaction parameter map pK : R → Rm1 × ... × 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 × nr matrix whose jth column is pK (r) T , where ρ(r) = j.The 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 T matrix is denoted by q.
We have the following useful Proposition: Proposition 7. Let (N , K) and (N , K ) be RIP-CFK systems.If T = T , then ψK = ψ K .
Proof.T = T ⇒ pK = p K for all K, K of the same type ⇒ IK = I K (by definition of interaction parameter map) ⇔ ψK = ψ 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 {K ∈ RIP K | ρ(r) = ρ(r ) ⇒ pK (r) = pK (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 ψK (x) = ψ K (x) ⇒ pK (r) = pK (r ) also holds, so that the corresponding sets are equal.Examples are PLK and P Y K 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 ρ(r) = ρ(r ) ⇒ pK (r) = pK (r ).For any other subset A of RIPK, we will denote {K ∈ RIP K | ρ(r) = ρ(r ) ⇒ pK (r) = pK (r )} with A-RDP (kinetics with reactantdetermined parameter maps).This notation is consistent with earlier ones since the corresponding letters there indicate the specific parameter maps, too.

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.

Key components of the MILP algorithm
The algorithm considers two CF systems: the original system (N , K) with N = (S , C , R) and the target system (N , K ) with N = (S , C , R ).The algorithm determines the corresponding network structure of the target system satisfying the linear conjugacy property.The two networks N and 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 R and K be known while R and K are to be obtained.The following are needed to be ascertained prior to the MILP implementation: where A k is the Laplacian map; • parameter > 0, that is set to be sufficiently small; and • parameter uij > 0, where i, j = 1, ..., m, i = j.Remark 8.Note that and u are introduced to ensure the correct structure of the linearly conjugate realization.(LC) 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 N has the correct structure.A dense linearly conjugate network can also be determined by considering the maximization problem analog.

MILP algorithm to CF systems
The optimal solution (if it exists) of the MILP would yield the matrix A b with the same structure as N , and the conjugacy constant vector c.The Laplacian map A k of the target network is computed as: where D = diag(e) and ej = c F •j , if complex j is a reactant of some reaction k 1, otherwise

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.

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.

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 , with I1 a PLK interaction map with kinetic order matrix F and (I2(x))j = m • dj • ej(x).Furthermore, the dissociation vectors dj (s.Definition 22) were organized in an r × 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.

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 X5 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:

Definition 4 .
The reaction vectors of a CRN (S , C , R) are the members of the set {y − y ∈ R S | (y, y ) ∈ R}.The stoichiometric subspace S of the CRN is the linear subspace of R S defined by S : span{y − y ∈ R S | (y, y ) ∈ R}.

Definition 18 .
If NR(y) is the number of CF-subsets of y, then 1 ≤ NR(y) ≤| ρ −1 (y) |.The reactant complex is a CF-node if NR(y) = 1, and an NF-node otherwise.It is a maximally NF-node if NR(y) =| ρ −1 (y) |> 1. Definition 19.The number N R of CF subsets of a CRN is the sum of NR(y) over all reactant complexes.

2 :INPUT1: reaction set with its kinetic values 3 : 4 : 5 :
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 ≤ (NR − nr) + l.Algorithm 1 CF-RM for RIP-NFK 1: procedure INITIAL OUTPUT1: reactant set, denote this by ρ(R) OUTPUT2: matrix ρ for the reactant map of the network (from OUTPUT1) procedure Identification of Branching Complexes 6: return complex is an NDK 22: procedure Generate RDK subsets for every NDK complexes 23: INPUT4: for every NDK node z (from OUTPUT4) 24:

Figure 2 .
Figure 2.: The "link-breaking" effect of CF-RM in the given CRN.

Figure 3 .
Figure 3.: The CRN after applying CF-RM to Schimtz's carbon cycle model in Figure 1.

Figure 5 .
Figure 5.: The graphs of the trajectories for (N * , K * ) and ( Ñ , K). M i represents a trajectory in the sparse realization.
The MILP algorithm finds a sparse linearly conjugate realization of the original network N .A sparse realization contains the minimum number of reactions, hence the associated objective function of the MILP model isMinimize m i,j=1 δij.(6.1)There are two sets of constraints in the model which indicate the linear conjugacy condition and desired structure of the network.

Figure 7 .
Figure 7.: The graphs of the trajectories for the original hill-type model and a linearly conjugate of sparse (first row) and dense (second row) realization.X i represents a trajectory in the sparse realization (shown in first row) and dense realization(shown in second row).

Table 1 .
1 lines : Key network number of a CF-RM transform.

Table 2
presents the network numbers of the N * .

Table 2 .
: Key network number of CRN N * of Schimtz's carbon cycle model.

Table 4 .
: List of variables used in the MILP Notation Description δ ij , i, j = 1, ..., m binary variable that keeps track of the presence of the reaction in the target netw [A b ] ij , i, j = 1, 2, ..., m kinetic matrix with the same structure as the target network n 21 )(k23 + (c3X3) n 23 ) −

Table 6 .
: List of symbols