Dynamical Equivalence and Linear Conjugacy of Chemical Reaction Networks: New Results and Methods

Dynamical Equivalence and Linear Conjugacy of Chemical Reaction Networks: New Results and Methods

Abstract

In the first part of this paper, we propose new optimization-based methods for the computation of preferred (dense, sparse, reversible, detailed and complex balanced) linearly conjugate reaction network structures with mass action dynamics. The developed methods are extensions of previously published results on dynamically equivalent reaction networks and are based on mixed-integer linear programming. As related theoretical contributions we show that (i) dense linearly conjugate networks define a unique super-structure for any positive diagonal state transformation if the set of chemical complexes is given, and (ii) the existence of linearly conjugate detailed balanced and complex balanced networks do not depend on the selection of equilibrium points. In the second part of the paper it is shown that determining dynamically equivalent realizations to a network that is structurally fixed but parametrically not can also be written and solved as a mixed-integer linear programming problem. Several examples illustrate the presented computation methods.


Keywords: chemical kinetics; stability theory; weak reversibility; linear programming; dynamical equivalence
AMS Subject Classifications: 80A30, 90C35.


1 Introduction

The mathematical study of chemical reaction networks is a rapidly growing field that has been applied recently to research problems in industrial chemistry, systems biology, gene regulation, and general nonlinear systems theory, among others [10, 23, 33]. There has also been significant theoretical work in the literature on such topics as persistence [1, 2, 3, 26, 32], multistability [8, 9, 31], monotonicity [6, 5], the global attractor conjecture for complex balanced systems [1, 2, 7, 12], lumping [15, 28, 39, 40], and conjugacy of reaction networks [11, 25].

One line of research which has been attracting increased attention has been that of determining when two reaction networks exhibit the same qualitative dynamics despite disparate network structure. In [11] and [34], the authors complete the question of what network structures can give rise to the same set of governing differential equations and therefore exhibit identical dynamics. This work was extended in [25] to networks which do not necessarily have the same set of differential equations but rather have trajectories related by a non-trivial linear transformation. Similar ground has also been touched in the papers [28] and [15] which deal with properties of linear lumpings, linear mappings which potentially reduce the dimension of the species set.

In this paper we look at the problem of algorithmically determining when a network is linearly conjugate to another network satisfying specified conditions. This problem was first addressed in [35] where the author presents a mixed-integer linear programming (MILP) algorithm capable of determining sparse and dense realizations, i.e. networks with the fewest and greatest number of reactions capable of generating the given dynamics. The algorithm was extended to complex and detailed balanced networks in [36], reversible networks in [37], weakly reversible networks in [38], and linear conjugate networks in [27], where a computationally more efficient procedure for determining weak reversibility was also presented.

We will continue in this paper the development and application of this MILP framework to linearly conjugate networks. In particular, we will give conditions which can be used to find sparse and dense linearly conjugate networks which are detailed and complex balanced, fully reversible, and which contain the greatest and fewest number of complexes. We will also expand the original MILP algorithm for finding sparse and dense realizations to find alternative realizations to a given reaction network when the network structure is fixed but the parameter values are not. Since this algorithm works without having to have the rate constants specified beforehand, this will allow us to answer questions about the reaction mechanism itself.

2 Background

In this section we present the terminology and notation relevant to chemical reaction networks and the main results from the literature upon which we will be building.

2.1 Chemical Reaction Networks

We will consider the chemical species or reactants of a network to be given by the set . The combined elements on the left-hand and right-hand side of a reaction are given by linear combinations of these species. These combined terms are called complexes and will be denoted by the set where

and the are nonnegative integers called the stoichiometric coefficients. We define the reaction set to be where the property will more commonly be denoted . To each we will associate a positive rate constant and to each we will set . The triplet will be called the chemical reaction network.

The above formulation naturally gives rise to a directed graph where the set of vertices is given by , the set of directed edges is given by , and the rate constants correspond to the weights of the edges from to . In the literature this has been termed the reaction graph of the network [24]. Since complexes may be involved in more than one reaction, as a product or a reactant, there is further graph theory we may consider. A linkage class is a maximally connected set of complexes, that is to say, two complexes are in the same linkage class if and only if there is a sequence of reactions in the reaction graph (of either direction) which connects them. A reaction network is called reversible if for any implies . A reaction network is called weakly reversible if for any implies there is some sequence of complexes such that .

A directed graph is called strongly connected if there exists a directed path from each vertex to every other vertex. A strongly connected component of a directed graph is a maximal set of vertices for which paths exists from each vertex in the set to every other vertex in the set. For a weakly reversible network, the linkage classes clearly correspond to the strongly connected components of the reaction graph.

Assuming mass-action kinetics, the dynamics of the specie concentrations over time is governed by the set of differential equations

(1)

where is the vector of reactant concentrations. The stoichiometric matrix contains entries and the Kirchhoff or kinetics matrix is given by

(2)

When we speak of the structure of a kinetics matrix, we will be referring to the distribution of positive and zero entries, which determines the network structure of the corresponding reaction graph. Finally, the mass-action vector is given by

(3)

2.2 Linearly Conjugate Networks

Under the assumption of mass-action kinetics, it is possible for the trajectories of two reaction networks and to be related by a linear transformation and therefore share many of the same qualitative features (e.g. number and stability of equilibria, persistence/extinction of species, dimensions of invariant spaces, etc.). This phenomenon was termed linear conjugacy in [25].

For completeness, we include the formal definition of linear conjugacy as presented in [25]. We will let denote the flow of (1) associated with and denote the flow of (1) associated with .

Definition 2.1.

We will say two chemical reaction networks and are linearly conjugate if there exists a linear function such that for all .

It is known that linear transformations can consist of at most positive scaling and reindexing of coordinates (Lemma 3.1, [25]). Linear conjugacy has been subsequently studied from a computational point of view in [27].

Linear conjugacy is a generalization of the concept of dynamical equivalence whereby two reaction networks with different topological network structure can generate the same exact set of differential equations (1). Two dynamically equivalent networks and are said to be alternate realizations of the kinetics (1), although it is sometimes preferable to say that is an alternative realization of or vice-versa. Since the case of two networks being realizations of the same kinetics is encompassed as a special case of linear conjugacy taking the transformation to be the identity, we will focus on linearly conjugate networks.

In general practice, we are given a network and asked to determine its dynamical behaviour. This is often a challenging problem; however, we may notice that the network behaves like one from a well-studied class of networks and therefore suspect a relationship which preserves key qualitative aspects of the dynamics. The theory of linear conjugacy can provide a powerful tool in analyzing such networks. If the network can be shown to be linearly conjugate to a network from the class of networks with understood dynamics, the dynamics of are transferred to .

This raises the question of how to find a linearly conjugate network when only the original network is given. This was studied in [27] where the authors built upon the linear programming algorithm introduced in [35]. We can impose that a network be linearly conjugate to our given network with the set of linear constraints

(4)

where , and the matrices and are given by:

(5)
(6)

The kinetics matrix for the network can by constructed from and by the relation

(7)

Finding a network satisfying (4) and then solving (7) is sufficient to determine a network which is linearly conjugate to via the transformation .

2.3 Sparse and Dense Linearly Conjugate Networks

In order to place the problem within a linear programming framework, we need to choose an objective function to optimize. An appropriate choice of such a function is not obvious and may vary depending on the intended application.

One particularly intuitive choice, which was introduced in [35] and has been widely used since, is to search for networks with the fewest and greatest number of reactions (sparse and dense networks, respectively). A sparse (respectively, dense) linearly conjugate network is given by a matrix satisfying (4) with the most (respectively, least) off-diagonal entries which are zeroes. Since the structure of and are the same, a correspondence between the non-zero off-diagonal entries in and a positive integer value can be made by considering the binary variables which will keep track of whether a reaction is ‘on’ or ‘off’, i.e. we have

where is sufficiently small and can be chosen the same as in (4), and where the symbol ’’ denotes the logical relation ‘if and only if’. These proposition logic constraints for the structure of a network can then be formulated as the following linear mixed-integer constraints (see, for example, [30]):

(8)

where for are appropriate upper bounds for the reaction rate coefficients. The number of reactions present in the network corresponding to is then given by the sum of the ’s so that the problem of determining a sparse network corresponds to satisfying the objective function

(9)

over the constraint sets (4) and (8). Finding a dense network corresponds to maximizing the same function, which can also be stated as a minimization problem as

(10)

It is known that, for trivial linear conjugacies, the structure of the dense realization contains the structure of all other trivial linear conjugacies as a subset (Theorem 3.1, [37]). We now prove the comparable result for non-trivial linear conjugacies.

Theorem 2.1.

Let be a chemical reaction network. Suppose that the reaction network is linearly conjugate to and dense. Suppose that is also linearly conjugate to . Then the directed unweighted graph of is a subset of the directed unweighted graph of .

Proof.

Assume and are linearly conjugate to , is dense in the space of networks which are linearly conjugate to , and contains a reaction not contained in .

Since both and are linearly conjugate to , we have by (4) that

and

where , and correspond to , and and correspond to .

Now consider . We have

where . Consequently, we have

(11)

On the other hand, we have

(12)

where . Combining (11) and (12) gives

where .

Since , and are diagonal matrices with positive entries on the diagonal, so is . This means that the network corresponding to is linearly conjugate to . We can readily see, however, that contains all of the reactions in both and . If contains a reaction not contained in then clearly has more reactions than which contradicts the assumption that is dense in the space of networks which are linearly conjugate to . The result follows. ∎

The following result follows immediately.

Corollary 2.1.

Let be a chemical reaction network. Then the structure of the unweighted directly graph of the dense reaction network which is linearly conjugate to is unique.

Proof.

This follows directly from Theorem 2.1. ∎

3 Computational Extensions of Linearly Conjugate Networks

In this section we will extend the optimization framework introduced in Section 2.3 to include complex balanced, reversible and detailed balanced networks, and to search for networks with the greatest and fewest number of complexes.

3.1 Weakly Reversible Networks

Weakly reversible networks are a particular important class of reaction networks because strong properties are known about their dynamics and equilibrium concentrations. Under a supplemental condition on the rate constants, weakly reversible networks are known to have complex balanced equilibrium concentrations and therefore exhibit all of the dynamical properties normally reserved for these networks [17, 22] (see Section 3.3 for further discussion of complex balanced networks).

Consequently, they are a primary candidate for the type of network we would like to find. The problem of determining if and when a chemical reaction network is linearly conjugate to a weakly reversible network was first considered in [38] and further refined in [27]. For convenience, we briefly recall the constraints published in [27] that guarantee the weak reversibility of the linearly conjugate network :

(13)

where is an auxiliary Kirchhoff matrix with the same structure as with appropriately scaled columns such that its kernel contains the -dimensional vector of all ones. In order to guarantee that the matrix has the same structure as and we also require that

(14)

3.2 Reversible Networks

In [36], an algorithm was presented which was capable of determining reversible reaction networks which are trivially linearly conjugate to a given reaction network. In this section, we present a simplified methodology and apply it to non-trivial linear conjugacies.

We recall that a network is reversible if for any implies . For the network , this is equivalent to the condition

for some sufficient small . This is in turn equivalent to

where , , as in Section 2.3. It follows that we can restrict our search space to reversible networks with the constraint set

(15)

A sparse or dense reversible network which is linearly conjugate to can be found by optimizing (9) or (10), respectively, over the constraint sets (4), (8) and (15).

3.3 Complex Balanced Systems

A particularly important class of chemical reaction networks are the complex balanced networks introduced in [24].

Definition 3.1.

An equilibrium concentration of the chemical reaction network is a complex balanced equilibrium concentration if

(16)

The network is called complex balanced if every equilibrium concentration is a complex balanced equilibrium concentration.

Many strong properties are known about complex balanced networks. In particular, it is known that complex balanced networks permit exactly one positive equilibrium concentration in each invariant space of the network and that this equilibrium concentration is locally asymptotically stable (Lemma 4C and Theorem 6A, [24]). Complex balanced systems are also known to be weakly reversible so that they are a subset of the weakly reversible networks considered in Section 3.1 (Theorem 2B, [22]).

The following result shows complex balancing is a system property depending on the structure and parameters of the network and not on the chosen equilibrium concentration.

Theorem 3.1 (Theorem 6A, [24]).

If a chemical reaction network is complex balanced at an equilibrium concentration then it is complex balanced at all of its equilibrium concentrations.

It should be noted that the complex balancing of a network is still dependent on the choice of rate constants. It is possible for a reaction network to be complex balanced for some choices of rate constants and not for others.

In [36], an algorithm was presented which was capable of determining sparse and dense complex balanced networks which are trivially linearly conjugate to a given network . This method required determining an equilibrium value of the network and then imposing the condition

on in accordance with (16). In this section, we extend these results to include non-trivial linearly conjugate networks.

Suppose that and are linearly conjugate via the transformation . In order to guarantee the network is complex balanced, according to (16) we require that

(17)

Since the equilibrium concentrations of and are related by the transformation , we have that the left-hand side of (17) can be rewritten

where we have made use of the form of the kinetics matrix of according to (7). The condition for complex balancing of the linearly conjugate network is therefore

(18)

where is as in (5). A sparse or dense complex balanced network which is linearly conjugate to can be found by optimizing (9) or (10), respectively, over the constraint sets (4), (8) and (18).

It should be noted that the optimization algorithm is less computationally exhausting than the corresponding algorithm for weak reversibility (13). This is because the matrix required in the general weakly reversible case is not required in the complex balancing condition; rather, it is sufficient to use the matrix . Consequently, there are fewer decision variables in the complex balancing algorithm. The pre-step that a be found satisfying may off-set this advantage, however, depending on the difficulty in solving .

It is unclear how the outcome of the algorithm depends on the choice of equilibrium concentration , which in general is not unique. The following result clarifies this dependence.

Theorem 3.2.

Suppose is linearly conjugate to with transformation matrix where and suppose is complex balanced at where and . Then is complex balanced at for all satisfying .

Proof.

Suppose trajectories of are related to trajectories of by the relationship where and .

Suppose that is complex balanced at where is an equilibrium concentration of . It follows from Theorem 3.1 that

(19)

Now consider an arbitrary equilibrium concentration of . Since and are linearly conjugate, it follows that . It follows that we have

It follows by the structure of that . From (19) we have that . In other words, is complex balanced at and we are done.

This result shows that when imposing the complex balancing constraint (18) on , it does not matter which equilibrium concentration of we choose. The feasible set of solutions (i.e. admissible networks) is the same.

3.4 Detailed Balanced Systems

In [36], the authors present an algorithm for determining detailed balanced networks which are trivially linearly conjugate to a given network. In this section we extend this algorithm to non-trivial linear conjugacies.

Definition 3.2.

An equilibrium concentration of the chemical reaction network is a detailed balanced equilibrium concentration if

(20)

The network is called detailed balanced if every equilibrium concentration is a complex balanced equilibrium concentration.

In other words, an equilibrium concentration is detailed balanced if the flow across each reaction is balanced by the flow across an opposite reaction at .

Suppose that and are linearly conjugate via the transformation . In order to guarantee the network is detailed balanced, according to (20) we require that

(21)

Since the equilibrium concentrations of and are related by the transformation , we have that

where we have made use of the form of the kinetics matrix of according to (7). The condition for detailed balancing of the linearly conjugate network is therefore

(22)

where is as in (5). A sparse or dense detailed balanced network which is linearly conjugate to can be found by optimizing (9) or (10), respectively, over the constraint sets (4), (8) and (22). Since the analogous result to Theorem 3.2 holds as a consequence of detailed balanced systems being a subset of complex balanced networks, we do not prove it here.

3.5 Minimal and Maximal Number of Complexes

We can also adapt to non-trivial linear conjugacies the algorithm introduced in [36] for determining a network with the fewest or greatest number of complexes from a fixed complex set which is trivially linearly conjugate to a given network .

In order to count the number of complexes in the network, we introduce the binary variables , , and consider the logical equations

(23)

for . In other words, takes on the value of one if and only if there is a reaction to or from the complex in the network; otherwise, it takes the value zero. For computational purposes, we reconsider (23) as

(24)

where . The linear constraints required to substantiate (24) are

(25)

We can now determine a network with the fewest or greatest number of complexes by optimizing the functions

(26)

or

(27)

respectively over the constraint sets (4) and (25). Further constraints can be imposed to restrict ourselves to specific classes of systems (e.g. complex balanced systems, reversible networks, etc.) although care has to be taken to ensure the structural constraints are still satisfied.

3.6 Examples

In this section, we present a few examples which illustrate the methodologies outlined so far.

Example 1: Consider the kinetic scheme

(28)

introduced in [27]. In that paper, it was shown that the kinetics could be generated by a reaction network involving the complex set

and that (28) has dynamics which is linearly conjugate to those generated by the sparse weakly reversible network given in Figure 1(a) (conjugacy constants ) and the dense weakly reversible network given in Figure 1(b) (conjugacy constants ).

Figure 1: Weakly reversible networks which are linearly conjugate to a network with the kinetics (28). The network in (a) is sparse while the networks in (b) and (c) are dense. The networks (a) and (c) are also complex balanced. (The parameter values in (b) and (c) have been rounded to three significant figures.)

The network in Figure 1(a) is complex balanced as a consequence of it being a zero deficiency weakly reversible network. It can be easily checked, however, that the network in Figure 1(b) is not complex balanced. We might wonder, therefore, what running the algorithm for a dense complex balanced network which is linearly conjugate to a network generating the kinetics (28) would produce.

Numerically, we can determine that an equilibrium concentration of (28) is . Running GLPK for a sparse network (9) over the constraints (4), (8) and (18) gives the network given in Figure 1(c) (conjugacy constants ). We notice that this network has the same structure as the weakly reversible network in Figure 1(b) and, furthermore, only differs in two rate constant values. It can be checked, however, that (c) is complex balanced while (b) is not.

Example 2: Consider the kinetic scheme

(29)

Using the indexing scheme introduced in [20] and more recently applied in the papers [36] and [27] we can construct a chemical reaction network capable of generating the dynamics (29) under the assumption of mass-action kinetics (1) which involves the complexes

It seems less than desirable, however, to consider a network of 18 complexes given the simplistic dynamics (29). We might wonder if there is a network with fewer complexes. Optimizing (26) over the constraint sets (4) and (25) gives the network

and the conjugacy constants .

In terms of understanding the qualitative dynamics of 29), this network is not particularly insightful. We might notice, however, that the net effect of all the reaction pathways leading out from is to create an and an at the expense of depleting . The complex , however, has not been considered in the procedure. Similarly, the reaction pathways leading out from generate an and an but the complex has not been included in the procedure. We might consider appending the procedure, therefore, to include and . Repeating the algorithm in GLPK gives the network

and the conjugacy constants . This is easily identified as the enzyme network

where an enzyme facilitates the transfer of a substrate into another substract and another enzyme facilitates the transfer back. This network was considered extensively in [3] and [4]. In particular, it was shown in [4] that for all rate constant values, the network possesses within each invariant space a unique positive equilibrium concentration which is globally asymptotically stable relative to that invariant space. It follows by the properties of linearly conjugate reaction networks that (29) inherits the same qualitative dynamics.

4 Structural Dynamical Equivalence

In this section we extend the computation procedure given in Section 3 for dynamical equivalence to the case of networks which are structurally fixed but have undetermined parameters.

4.1 Dynamical Equivalence

The MILP framework outlined so far requires that the rate constants be specified for the network . Consequently, when we search for networks which are linearly conjugate to a given network , we are really asking if there are networks which are linearly conjugate for a given choice of parameter values.

For networks where the dynamical behaviour is heavily dependent on the chosen rate constants, however, it is possible that certain behaviours are being overlooked by poor rate constant selection. There are networks, for instance, which are known to be linearly conjugate to weakly reversible networks or complex balanced networks for certain values of the rate constants but not for others (see Examples 2 and 3 of [25]). In such cases, if the rate constants are not carefully chosen, the algorithm may overlook these networks and we would not realize that the mechanism shares characteristics with these other networks.

Therefore, we now change the problem setup by fixing only the structure of the initial network but not the parameter values. We will show that this problem class can also be casted to the framework of MILP. This would remove the above mentioned limits of using a fully specified initial network model.

The conditions for dynamical equivalence, keeping the entries of both and general, are

(30)

Note that in (30) both the off-diagonal entries of and are now decision variables.

As in Section 2.3, we want to keep track of the structure of . The conditions corresponding to (8) for the matrix are

(31)

As before, the binary variables keep track of whether a reaction is in the network or not and thus are capable of counting the number of reactions in .

We also, however, want to permit to have a variable rate constant values within a fixed network structure. In order to fix this network structure, we introduce the binary variables , , , and the logical equations

(32)

for some . In other words, the ’s keep track of whether the reaction is in the network . The conditions required to allow the entries of to vary independently within this pre-defined structure are

(33)

where

(34)

Further constraints can be implemented to search through subspaces of the parameter spaces for alternative realizations. For instance, if we suspect the reaction rate for the reaction is slaved to that of we can add to the procedure, etc.

The conditions (30), (31) and (33) can be combined with the structural conditions for reversibility (15) and weak reversibility (13 and 14) and the objective functions (9) and (10) to search over the parameter space of for sparse and dense alternative realizations which satisfy these further structural constraints.

4.2 Complex Balanced Realizations

It is also desirable to explore the parameter space of for alternative realizations which are complex balanced. The linear constraints (22) and (18), however, cannot be used in the parameter-independent case since the required equilibrium concentrations depend on the rate constants for which are not specified.

We might expect, however, since all weakly reversible networks are complex balanced for some choice of rate constants, that if a network has a weakly reversible alternative realization for some other choice of rate constants then it also has a complex balanced alternative realization for some choice of rate constants. That is to say, in the parameter-independent optimization procedure weak reversibility is sufficient to demonstrate complex balancing. The main result of this subsection guarantees this (Theorem 4.3).

First, however, we need the following results about weakly reversible networks.

Theorem 4.1 (Theorem 3.1, [19]; Proposition 4.1, [16]).

Let be a kinetics matrix and let , denote the support of the linkage class. Then the reaction graph corresponding to is weakly reversible if and only if there is a basis of ker, , such that, for ,

Theorem 4.2 (Theorem 1, [13]).

Under the assumption of mass-action kinetics, weakly reversible chemical reaction networks possess at least one positive equilibrium concentration within each positive invariant space of the system.

An immediate consequence of Theorem 4.1 is that a network is weakly reversible if and only if there is a vector in the kernel of . We will exploit this fact in the next result.

Theorem 4.3.

Suppose there is a choice of rate constants such that the network is dynamically equivalent to and is weakly reversible. Then there exists a choice of rate constants such that the network is dynamically equivalent to where is complex balanced. Furthermore can be selected to have the same structure as .

Proof.

Let be the kinetics matrix associated with and be the kinetics matrix associated with , and suppose that is weakly reversible. Let denote the positive vector in ker guaranteed to exist by Theorem 4.1 and be any positive equilibrium concentration of (1) guaranteed to exist by Theorem 4.2.

We now define a new network with the associated kinetics matrix