Dynamical equivalence

From Mathematics of Reaction Networks
Jump to: navigation, search

Two mass action systems (\mathcal{S}, \mathcal{C}, \mathcal{R},\mathcal{K}) and (\mathcal{S},\mathcal{C}',\mathcal{R}',\mathcal{K}') are said to be dynamically equivalent if the differential equations produced under the assumption of mass action kinetics are identical. Two mass action systems which given rise to the same differential equation set are sometimes said to be alternative realizations of the same kinetics.

It has been known since at least the early 1970s that two mass action systems with distinct rate constants and network structure can given rise to the same dynamics. Frederick J. Krambeck considered such systems arising from detailed balancing[1] (which he called non-uniqueness of parameters) and Fritz Horn and Roy Jackson considered it for a class of weakly reversible networks[2] (which they called macro-equivalence). A thorough study of the capacity of networks to be dynamically equivalence was conducted by Gheorghe Craciun and Casian Pantea[3] (they called the phenomenon counfoundability of mass action systems).

Recent research on dynamical equivalence of mass action systems has been focused on two primary problems:

  • Identifiability of network structure - The reality that mass action systems with different underlying network structure may give rise to the same dynamics provides a limitation for the inverse problem of chemical kinetics. Specifically, even if we have perfect knowledge about the dynamics of a chemical kinetic system, we may not be able to uniquely determine the network structure which gave rise to the dynamics. This is a significant concern and limitation in application as we typically only have measurable information at the end of a process and must infer the network structure. Exploring the limits of permissible networks with dynamically equivalent kinetics has consequently been a focus in chemical reaction network theory[4][3][5][6].
  • Correspondence of dynamics - It is possible for two mass action systems to be dynamically equivalent even when only one of them has a typically desirable feature about either their dynamics (e.g. detailed balancing, complex balancing, etc.) or their underlying network structure (e.g. reversibility, weak reversibility, low deficiency, etc.). In these cases, the dynamical property guaranteed by the desirable systems applies to the system which does not have the property. In other words, the scope of mass action systems with the desirable property is expanded. Determining conditions which expand the scope of applicability of classical results (e.g. deficiency theory, detailed balancing, complex balancing, etc.) has been a focus of research [7][8] and in recent years specifically has taken a turn toward computational methods[9][10][11][12][13][14][15]

Identifiability of chemical reaction networks

Correspondence of dynamics

It is possible for two mass action systems (\mathcal{S}, \mathcal{C},\mathcal{R},\mathcal{K}) and (\mathcal{S},\mathcal{C}',\mathcal{R}',\mathcal{K}') to exhibit qualitative identical behavior despite difference in network structure, and even subtle differences in the underlying differential equations governing the system. It is particularly noteworthy that this may hold despite significant differences in the connectedness and reversibility of the reaction graphs of the chemical reaction networks underlying the systems. That is to say, one of the networks underlying the mass action systems involved may have a desirable structural property (weak reversibility, endotacticity, single linkage class, etc.) while the other may not.

The realization that two systems exhibit the same dynamics can be a powerful analysis tool when results are known which guarantee one of the systems has known dynamic properties (e.g. locally stable dynamics, multistationarity, persistence, etc.). These tools can therefore be used to expand the scope of existing theory to systems for which the underlying network would not otherwise guarantee such behavior.

Dynamical equivalence

It is possible for two mass action systems (\mathcal{S}, \mathcal{C},\mathcal{R},\mathcal{K}) and (\mathcal{S},\mathcal{C}',\mathcal{R}',\mathcal{K}') to give rise to identical systems systems of differential equations despite differences in the complex set \mathcal{C}, the reaction set \mathcal{R}, and/or the kinetic rate constants \mathcal{K}.

For example, consider the four-complex network considered by Fritz Horn and Roy Jackson[2]:

2\mathcal{A}_1 + \mathcal{A}_2 \; \stackrel{1}{\rightarrow} \; 3\mathcal{A}_1 \\
{}^{\epsilon} \uparrow \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \downarrow {}_{\epsilon} \\
3\mathcal{A}_2 \; \stackrel{1}{\leftarrow} \; \mathcal{A}_1 + 2\mathcal{A}_2

The network is weakly reversible but has a deficiency of two so that the mass action system does not fall within the scope of the Deficiency Zero Theorem. The conditions of the General Deficiency Theorem guarantee that the mass action system is complex balanced if and only if \epsilon = 1. Outside of this value, deficiency theory is silent.

Consider, however, the reaction vector associated with the reaction from the complex 2\mathcal{A}_1 + \mathcal{A}_2. For \epsilon > 1, we can rescale and split the mass action term corresponding to this reaction according to

\left[ \begin{array}{c} 1 \\ -1 \end{array} \right] c_1^2c_2 = \left( \frac{1+2\epsilon}{3} \left[ \begin{array}{c} 1 \\ -1 \end{array} \right] + \frac{\epsilon -1}{3} \left[ \begin{array}{c} -2 \\ 2 \end{array} \right] \right) c_1^2c_2.

In other words, the reaction can be split into two, one linking to 3\mathcal{A}_1 and one linking to 3\mathcal{A}_2. Performing the same analysis on the reaction from \mathcal{A}_1+2\mathcal{A}_2 yields the dynamically equivalent system

2\mathcal{A}_1 + \mathcal{A}_2 \; \stackrel{(1+2\epsilon)/3}{\longrightarrow} \; 3\mathcal{A}_1 \\
\frac{\epsilon-1}{3} \downarrow \uparrow \epsilon \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \epsilon \downarrow \uparrow \frac{\epsilon-1}{3} \\
3\mathcal{A}_2 \; \stackrel{(1+2\epsilon)/3}{\longleftarrow} \; \mathcal{A}_1 + 2\mathcal{A}_2

for which the equilibrium concentration \bar{c}_1=\bar{c}_2 is always a complex balanced equilibrium concentration. In other words, we can show that the original network behaves like a complex balanced system for all \epsilon > 1 by rewriting the system in a way that it corresponds to a system with different underlying network structure.

Consider similarly splitting the reaction from 3 \mathcal{A}_1. For 1/2 \leq \epsilon < 1 we have

\epsilon \left[ \begin{array}{c} -2 \\ 2 \end{array} \right]c_1^3 = \left( \left( 2\epsilon -1 \right) \left[ \begin{array}{c} -2 \\ 2 \end{array} \right] + \left( 2 - 2 \epsilon \right) \left[ \begin{array}{c} -1 \\ 1 \end{array} \right] \right) c_1^3.

In other words, we can split the reaction into two, one linking to \mathcal{A}_1 + 2\mathcal{A}_2 and one linking to 2 \mathcal{A}_1 + \mathcal{A}_2. Performing the analogous procedure on the reaction originating from 3\mathcal{A}_2 gives us the dynamically equivalent network

\begin{array}{c} 2 \mathcal{A}_1 + \mathcal{A}_2 \; \; \stackrel{2-2\epsilon}{\leftrightarrows} \; \; 3\mathcal{A}_1 \\ {}^{2\epsilon - 1} \uparrow \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \downarrow {}_{2\epsilon-1} \\ 3\mathcal{A}_2 \; \; \stackrel{2-2\epsilon}{\rightleftarrows} \; \; \mathcal{A}_1 + 2\mathcal{A}_2\end{array}

where the unlabelled reactions have rate constant 1. The equilibrium concentration \bar{c}_1 = \bar{c}_2 is always a complex balanced equilibrium concentration for this system. Piecing this together with the previous decomposition, we have that by restructuring the mass action equations, we can show that the network behaves like a complex balanced mass action system (that is to say, exhibits locally stable dynamics) for all \epsilon \geq 1/2 despite being only complex balanced for the value \epsilon = 1. This is a significant improvement over the results guaranteed by consider the original formulation of the network, although it does fall short of capturing the whole range of parameter values capable of producing locally stable dynamics, which is known to be \epsilon \geq 1/6.

Linear conjugacy

Two mass action systems (\mathcal{S},\mathcal{C},\mathcal{R},\mathcal{K}) and (\mathcal{S},\mathcal{C}',\mathcal{R}',\mathcal{K}') for which the flows of the corresponding differential equations are related by a non-trivial linear transformation are said to be linearly conjugate. Formally, suppose that the mass action system (\mathcal{S},\mathcal{C},\mathcal{R},\mathcal{K}) has flow \Phi(\vec{x}_0,t) and the mass action system corresponding to (\mathcal{S},\mathcal{C}',\mathcal{R}',\mathcal{K}') has flow \Psi(\vec{y}_0,t). The two systems are said to be linearly conjugate to one another if there exists a linear mapping T: \mathbb{R}^m_{\geq 0} \mapsto \mathbb{R}^m_{\geq 0} such that T \Phi(\vec{x}_0,t) = \Psi(T\vec{x}_0,t) for all t \geq 0.

Linearly conjugacy was introduced by Matthew Johnston and David Siegel in [16] as a generalization of dynamical equivalence. Dynamical equivalence can be thought of as a trivial linear conjugacy; that is to say, a linearly conjugate relationship where the conjugate map is the identity. Linear conjugacy is also related to the concept of lumping.

For example, consider the chemical reaction networks[15]

\mathcal{N}: \; \; \; \; \; 
\displaystyle{\mathcal{A}_1 \; \stackrel{k_1}{\longrightarrow} \; \mathcal{A}_2}\\
\displaystyle{2\mathcal{A}_2 \; \stackrel{k_2}{\longrightarrow} \; 2\mathcal{A}_1}
\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \mathcal{N}': \; \; \; \; \; 
\displaystyle{\mathcal{A}_1 \; \stackrel{\tilde{k}_1}{\longrightarrow} \; 2\mathcal{A}_2}\\
\displaystyle{2\mathcal{A}_2 \; \stackrel{\tilde{k}_2}{\longrightarrow} \; \mathcal{A}_1.}

Under the assumption of mass action kinetics, these networks give rise to the following governing set of ordinary differential equations

\displaystyle{\frac{dx_1}{dt}} & = -k_1 x_1 + 2k_2 x_2\\
\displaystyle{\frac{dx_2}{dt}} & = k_1 x_1 - 2k_2 x_2 
\; \; \; \; \; \; \; \; \; \; \mbox{and} \; \; \; \; \; \; \; \; \; \;
\displaystyle{\frac{dy_1}{dt}} & = -\tilde{k}_1 y_1 +\tilde{k}_2 y_2\\
\displaystyle{\frac{dy_2}{dt}} & = 2\tilde{k}_2 y_1 - 2 \tilde{k}_2 y_2.

It is clear that these networks are not dynamically equivalent since they do not have an identical set of governing differential equations. For any set of rate constants, however, trajectories for the system corresponding to \mathcal{N} can be related to that of \mathcal{N}' by the linear relationship y_1(t)=x_1(t), y_2(t)=2x_2(t). Since \mathcal{N}' is a reversible network with zero deficiency, it follows by the Deficiency zero theorem that \mathcal{N}' exhibits locally stable dynamics. It follows by the nature of the conjugacy mapping that \mathcal{N} exhibits the same qualitative dynamics.

It is known that the linear mapping T may preserve mass action kinetics if and only if it has precisely one nonzero entry in each row and column, and that entry is positive (Lemma 3.1 [16]). In other words, the mapping may only permute and rescale the species. It is also common to restrict permutation so that T is a diagonal matrix with positive weighted entries on its diagonal.

Computational approaches


  1. Frederick J. Krambeck, The mathematical structure of chemical kinetics in homogeneous single-phase systems, Arch. Rational Mech. Anal., 38:317--347, 1970
  2. 2.0 2.1 Fritz Horn and Roy Jackson, General mass action kinetics, Arch. Ration. Mech. Anal., 47:187--194, 1972
  3. 3.0 3.1 Gheorghe Craciun and Casian Pantea, Identifiability of chemical reaction networks, J. Math Chem., 44(1):244--259, 2008
  4. V. Hars and Janos Toth, On the inverse problem of reaction kinetics, Coll. Math. Soc. J. Bolyai, 30:363--379, 1981
  5. Gheorghe Craciun, Casian Pantea, and Grzegorz Rempala, Algebraic methods for inferring biochemical networks: a maximum likelihood approach, Compute. Biol. Chem., 33(5):361--367, 2009
  6. Gheorghe Craciun, Jaejik Kim, Casian Pantea, and Grzegorz A. Rempala, Statistical Model for Biochemical Networks Inference, to appear in Communications in Statistics: Simulation and Computation, 2012
  7. Einstein D. Averbukh, Some equivalent kinetic functions of macrodynamic equations, Automation and Remote Control, 55(11):1694--1698, 1994
  8. Einstein D. Averbukh, Complex balanced kinetic functions in inverse problems of chemical kinetics, Automation and Remote Control, 55(12):1723--1732, 1994
  9. Gabor Szederkenyi, Computing sparse and dense realizations of reaction kinetic systems, J. Math. Chem., 47:551--568, 2010
  10. Gabor Szederkenyi and Katalin Hangos, Finding complex balanced and detailed balanced realizations of chemical reaction networks, J. Math. Chem., 49:1163--1179, 2011
  11. Gabor Szederkenyi, Katalin Hangos, and Tamas Peni, Maximal and minimal realizations of chemical kinetics systems: computation and properties.MATCH Commun. Math. Comput. Chem., 65:309--332, 2011
  12. Gabor Szederkenyi, Katalin Hangos, and Zsolt Tuza, Finding weakly reversible realizations of chemical reaction networks using optimization. MATCH Commun. Math. Comput. Chem., 67:193--212, 2012
  13. Matthew D. Johnston, David Siegel, and Gabor Szederkenyi, A linear programming approach to weak reversibility and linear conjugacy of chemical reaction networks, J. Math. Chem., 50(1):274--288, 2012
  14. Matthew D. Johnston, David Siegel, and Gabor Szederkenyi, Dynamical equivalence and linear conjugacy of chemical reaction networks: New results and methods, to appear in MATCH Commun. Math. Comput. Chem., 2012
  15. 15.0 15.1 Matthew D. Johnston, David Siegel, and Gabor Szederkenyi, Computing weakly reversible linearly conjugate chemical reaction networks with minimal deficiency, to appear in Math. Biosci., 2012
  16. 16.0 16.1 Matthew D. Johnston and David Siegel, Linear conjugacy of chemical reaction networks, J. Math. Chem., 49(7):1263--1282, 2011.