A computational approach to the structural analysis of uncertain kinetic systems

TL;DR

This paper introduces a computational framework for analyzing uncertain kinetic systems, focusing on reaction graph structures within models where coefficients vary within a polytopic set, enabling efficient structural analysis.

Contribution

It presents a polynomial-time method to compute the dense realization and all possible reaction graph structures of uncertain kinetic models.

Findings

Dense realization contains the maximum number of reactions.

Core reactions are present in all reaction graphs.

Algorithm for enumerating all reaction graph structures.

Abstract

A computation-oriented representation of uncertain kinetic systems is introduced and analysed in this paper. It is assumed that the monomial coefficients of the ODEs belong to a polytopic set, which defines a set of dynamical systems for an uncertain model. An optimization-based computation model is proposed for the structural analysis of uncertain models. It is shown that the so-called dense realization containing the maximum number of reactions (directed edges) is computable in polynomial time, and it forms a super-structure among all the possible reaction graphs corresponding to an uncertain kinetic model, assuming a fixed set of complexes. The set of core reactions present in all reaction graphs of an uncertain model is also studied. Most importantly, an algorithm is proposed to compute all possible reaction graph structures for an uncertain kinetic model.

Figures

Equations

$$\dot{x}=M\cdot \phi (x)$$

$$[Y{]}_{ij}={\alpha }_{ji}i\in \{1,\dots ,n\},j\in \{1,\dots ,m\}$$

$$[A_k{]}_{ij}=\{\begin{array}{ll}{k}_{ji}& \text{ if }i\ne j\\ -\sum _{l=1,l\ne i}^m{k}_{il}& \text{ if }i=j\end{array}i,j\in \{1,\dots ,m\}$$

$$\dot{x}=Y\cdot A_k\cdot {\psi }^Y(x)$$

$${\psi }_j^Y(x)=\prod _{i=1}^n{x}_i^{[Y{]}_{ij}},j\in \{1,\dots ,m\}$$

$$M\cdot \phi (x)=Y\cdot A_k\cdot {\psi }^Y(x)$$

$$M=Y\cdot A_k$$

$$\mathcal{P}=\{\tilde{M}\in {\mathbb{R}}^{nm}|{\tilde{M}}^{\top }\cdot n_i\le b_i,1\le i\le q\}\}$$

$$[Y{]}_{ij}=0⟹[M{]}_{ij}\ge 0,i\in \{1,\dots ,n\},j\in \{1,\dots ,m\}$$

$${\tilde{M}}^{\top }\cdot {\alpha }_i+{\widetilde{A_k}}^{\top }\cdot {\beta }_i\le d_i$$

$${\tilde{M}}^{\top }\cdot n_i\le b_i,$$

$$M=Y\cdot A_k,$$

$$[A_k{]}_{ij}\ge 0,$$

$$\sum _{j=1}^m[A_k{]}_{ij}=0,$$

$$T=c\cdot D+(1-c)\cdot Rc\in (0,1)$$

$$b_i=\{\begin{array}{ll}1& \text{ if }i\in B\\ 0& \text{ if }i\notin B\end{array}$$

$${\dot{x}}_1=3c_1\cdot x_2^3-c_2\cdot x_1^3$$

$${\dot{x}}_2=-3c_1\cdot x_2^3+c_2\cdot x_1^3,$$

$${\widetilde{M_{}}}^{\top }\cdot e_l^6\le (1+{\gamma }_l)\cdot [M{]}_{ij}$$

$${\widetilde{M_{}}}^{\top }\cdot (-e_l^6)\le ({\rho }_l-1)\cdot [M{]}_{ij}$$

$$Y\cdot \left[\begin{array}{c}-k_{12}-k_{13}\\ k_{12}\\ k_{13}\end{array}\right]=\left[\begin{array}{c}3c_1\\ -3c_1\end{array}\right]k_{12},k_{13}\in {\mathbb{R}}^{+}⟹k_{12}\in [0,c_1],k_{13}=\frac{3}{2}{c}_1-\frac{3}{2}{k}_{12}$$

$$Y\cdot \left[\begin{array}{c}{k}_{21}\\ -k_{21}-k_{23}\\ k_{23}\end{array}\right]=\left[\begin{array}{c}-c_2\\ c_2\end{array}\right]k_{21},k_{23}\in {\mathbb{R}}^{+}⟹k_{21}\in (0,\frac{c_2}{3}),k_{23}=c_2-3k_{21}$$

$$Y\cdot \left[\begin{array}{c}{k}_{31}\\ k_{32}\\ -k_{31}-k_{32}\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]k_{31},k_{32}\in {\mathbb{R}}^{+}⟹k_{31}\in {\mathbb{R}}^{+},k_{32}=2k_{31}$$

$${\widetilde{M_{}}}^{\top }\cdot (-e_1^6)$$

$${\widetilde{M_{}}}^{\top }\cdot (-e_5^6)$$

$${\widetilde{M_{}}}^{\top }\cdot e_3^6$$

$${\widetilde{M_{}}}^{\top }\cdot e_6^6$$

$${\widetilde{M_{}}}^{\top }\cdot [1,1,0,1,1,0{]}^{\top }$$

$${\widetilde{M_{}}}^{\top }\cdot [0,-1,0,-1,0,0{]}^{\top }$$

$${\widetilde{M_{}}}^{\top }\cdot [-1,0,0,0,1,0{]}^{\top }$$

$$Y$$

$$A_k$$

Full text

Keywords: reaction networks, uncertain models, reaction graphs, algorithms, convex optimization

1 Introduction

Kinetic models in the form of nonlinear ordinary differential equations are widely used for describing time-varying physico-chemical quantities in (bio-)chemical environments [46]. Moreover, the kinetic system class is dynamically rich enough to characterize general nonlinear behaviour in other application fields as well, particularly where the state variables are nonnegative and the model has a networked structure, such as in the modelling of process systems, population or disease dynamics, or even transportation processes [14, 5, 38]. In biochemical applications, the exact values (or even sharp estimates) of the model parameters are often not known, making the models uncertain [6]. Even when we have measurements of sufficient quantity and quality, the lack of structural or practical identifiability may result in highly uncertain models even with the most sophisticated estimation methods [10, 34, 8, 7]. This inherent uncertainty was a key factor in the development of Chemical Reaction Network Theory (CRNT), where (among other goals) a primary interest is to study the relations between the network structure and the qualitative properties of the corresponding dynamics, preferably without the precise knowledge of model parameters. From the earlier results of CRNT, we have to mention the well-known Deficiency One and Deficiency Zero Theorems [12, 32] opening the way towards a structure-based (essentially parameter-free) dynamical analysis of biological networks. Recent particularly important findings in this area are the identification of biologically plausible structural sources of absolute concentration robustness [31], and the proof of the Global Attractor Conjecture [1, 9].

The efficient treatment of uncertain quantitative models is a fundamental task in mathematics, physics, (bio)chemistry and in related engineering fields [4, 13]. An important early result is [15], where the solutions of linear compartmental systems are studied with uncertain flow rates that are assumed to belong to known intervals. In [22] a probabilistic framework is proposed for the representation and analysis of uncertain kinetic systems. In [28] an analytical expression is computed for the temperature dependence of the uncertainty of reaction rate coefficients, and a method is proposed for computing the covariance matrix and the joint probability density function of the Arrhenius parameters. A recent outstanding result is [29], where a deterministic computation interpolation scheme for uncertain reaction network models is proposed, which is able to handle large-scale models with hundreds of species and kinetic parameters.

The description of model uncertainties using convex sets is often a computationally appealing way of solving model analysis, estimation or control problems [40, 3]. From the numerous applications, we mention here only a few selected works from different fields. In [41], a stabilization scheme was given for nonlinear control system models, where the uncertain coefficients of smooth basis functions in the system equations are assumed to form a polytopic set. An interval representation of fluxes in metabolic networks was introduced in [26], which enables the computation of the $\alpha$-spectrum even from an uncertain flux distribution. In [24], a nonlinear feedback design method is proposed which is able to robustly stabilize parametrically uncertain kinetic systems using the convexity of the constraint ensuring the complex balance property. Recently, a new approach was given for the stability analysis of general Lotka-Volterra models with polytopic parameter uncertainties in [2].

It is known from the fundamental dogma of chemical kinetics that the reaction graph structure corresponding to a kinetic ODE-model is generally non-unique, even in the case when the rate coefficients are assumed to be known [16, 46, 30]. This property is usually called dynamical equivalence, macro-equivalence or confoundability in the literature [10, 16]. The first solution to the inverse problem, namely the construction of one possible reaction network (called the canonical network) for a given set of kinetic differential equations was described in [17]. The notion of dynamical equivalence was extended by introducing linear conjugacy of kinetic systems in [18] allowing a positive diagonal transformation between the solutions of the kinetic differential equations. The simple factorization of kinetic models containing the Laplacian matrix of the reaction graph allows the development of efficient methods in various optimization frameworks for computing reaction networks realizing or linearly conjugate to a given dynamics with preferred properties such as density/sparsity [33], weak reversibility [20], complex or detailed balance [35], minimal or zero deficiency [21, 23]. Using the superstructure property of the so-called dense realizations, it is possible to algorithmically generate all possible reaction graph structures corresponding to linearly conjugate realizations of a kinetic dynamics [45, 43].

Even if the monomials of a kinetic system are known, the parameters (i.e., the monomial coefficients) are often uncertain in practice. For example, one may consider the situation when a kinetic polynomial ODE model with fixed structure is identified from noisy measurement data. In such a case, using the covariance matrix of the estimates and the nonnegativity/kinetic constraints for the system model, we can define a simple interval-based (see, e.g. [26]), or more general (e.g., polytopic or ellipsoidal) uncertain model [25, 37]. Based on the above, the goal of this paper is to extend and illustrate previously introduced notions, computational models and algorithms for kinetic systems with polytopic uncertainty.

2 Notations and computational background

In this section we summarize the basic notions of kinetic polynomial systems and the generalized model defined with uncertain parameters.

The applied general notations are listed below:

[TABLE]

2.1 Kinetic polynomial systems and their models

Nonnegative polynomial systems are defined in the following general form:

[TABLE]

where $x:\mathbb{R}\rightarrow\mathbb{R}_{+}^{n}$ is a nonnegative valued function, $M\in\mathbb{R}^{n\times p}$ is a coefficient matrix and $\varphi:\mathbb{R}_{+}^{n}\rightarrow\mathbb{R}_{+}^{p}$ is a monomial-type vector-mapping. The invariance of the nonnegative orthant with respect to the dynamics (1) can be ensured by prescribing sign conditions for the entries of matrix $M$ depending on the exponents of $\varphi$, see [46, 14].

In this paper, we treat kinetic models as a general nonlinear system class that is suitable for the description of biochemical reaction networks. Hence, we do not require that all models belonging to the studied class are actually chemically realizable. Several physically or chemically relevant properties such as component mass conservation, detailed or complex balance can be ensured by adding further constraints to the computations (see, e.g. [19]).

Definition 2.1**.**

*A chemical reaction network (CRN) can be characterized by three sets [11, 12].

species: $\mathcal{S}=\{X_{i}\ |\ i\in\{1,\ldots,n\}\}$

complexes: $\mathcal{C}=\{C_{j}=\sum\limits_{i=1}^{n}\alpha_{ji}X_{i}\ |\ \alpha_{ji}\in\mathbb{N},\ j\in\{1,\ldots,m\},\ i\in\{1,\ldots,n\}\}$

reactions: $\mathcal{R}\subseteq\{(C_{i},C_{j})\ |\ C_{i},C_{j}\in\mathcal{C}\}$

For all $i,j\in\{1,\ldots m\}$, $i\neq j$ the reaction $C_{i}\rightarrow C_{j}$ is represented by the ordered pair $(C_{i},C_{j})$, and it is described by a nonnegative real number $k_{ij}\in\mathbb{R}^{+}$ called reaction rate coefficient. The reaction $C_{i}\rightarrow C_{j}$ is present in the reaction network if and only if $k_{ij}$ is strictly positive.*

The relation between species and complexes is described by the complex composition matrix $Y\in\mathbb{R}^{n\times m}$, the columns of which correspond to the complexes, i.e.

[TABLE]

The presence of the reactions in the CRN is defined through the rate coefficients as the off-diagonal entries of the Kirchhoff matrix $A_{k}\in\mathbb{R}^{m\times m}$ which is a Metzler compartmental matrix with zero column-sums. Its entries are defined as:

[TABLE]

According to this notation, the reaction $C_{i}\rightarrow C_{j}$ takes place in the reaction network if and only if $[A_{k}]_{ji}$ is positive, and $[A_{k}]_{ji}=0$ implies that $(C_{i},C_{j})\notin\mathcal{R}$. Since a chemical reaction network is uniquely characterized by the matrices $Y$ and $A_{k}$, we refer to a CRN by the corresponding pair $(Y,A_{k})$.

If mass action kinetics is assumed, the equations governing the dynamics of the concentrations of the species in the CRN defined by the function $x:\mathbb{R}\rightarrow\mathbb{R}^{n}_{+}$ can be written in the form:

[TABLE]

where $\psi^{Y}:\mathbb{R}_{+}^{n}\rightarrow\mathbb{R}_{+}^{m}$ is the monomial function of the CRN with coordinate functions

[TABLE]

The nonnegative polynomial system (1) is called a kinetic system if there exists a reaction network $(Y,A_{k})$ so that its dynamics satisfies the equation [46]:

[TABLE]

As it has been mentioned in the Introduction, reaction networks with different sets of complexes and reactions may be governed by the same dynamics. If Equation (6) is fulfilled, then the CRN $(Y,A_{k})$ is called a dynamically equivalent realization of the kinetic system (1).

The description of the polynomial system (1) can be transformed so that the monomial function $\varphi$ is equal to $\psi^{Y}$ (and $p=m$ holds) while the described dynamics remains the same. After the transformation and simplification based on the properties of polynomials, Equation (6) can be simplified to:

[TABLE]

Reaction networks have another representation, which is more suitable for illustrating the structural properties. It is a weighted directed graph $G(V,E)$ called the Feinberg-Horn-Jackson graph or reaction graph for brevity [46]. The complexes are represented by the vertices, and the reactions by the edges. Let the vertices $v_{i}$ and $v_{j}$ correspond to the complexes $C_{i}$ and $C_{j}$, respectively. Then there is a directed edge $v_{i}v_{j}\in V(G)$ with weight $k_{ij}$ if and only if the reaction $C_{i}\rightarrow C_{j}$ takes place in the CRN.

2.2 Uncertain kinetic systems

For the uncertainty modelling, we assume that the monomial coefficients in matrix $M$ are constant but uncertain, and they belong to an $n\cdot m$ dimensional polyhedron.

Remark 2.2**.**

In previous sections the set of uncertain parameters is noted as a polytope or a polytopic set, but from now on we use the notion of a polyhedron as well. The former one is defined as the convex hull of its vertices, while the latter one is the intersection of halfspaces, and the two definitions are not equivalent in general. However, in the examined problems it is assumed that the parameters of the kinetic models are bounded, and a bounded polyhedron is equivalent to a bounded polytope.

We represent the matrix $M$ as a point denoted by $\widetilde{M}$ in the Euclidean space $\mathbb{R}^{nm}$. In the uncertain model it is assumed that the possible points $\widetilde{M}$ are all the points of a closed convex polyhedron $\mathcal{P}$, which is defined as the intersection of $q$ halfspaces. The boundaries of the halfspaces are hyperplanes with normal vectors $n_{1},\ldots,n_{q}\in\mathbb{R}^{nm}$ and constants $b_{1},\ldots,b_{q}\in\mathbb{R}$. Applying these notations, the polyhedron $\mathcal{P}$ can be described by a linear inequality system as

[TABLE]

For the characterization of the polyhedron $\mathcal{P}$ not only the possible values of the parameters should be considered, but also the kinetic property of the polynomial system. This can be ensured (see [17]) by prescribing the sign pattern of the matrix $M$ as follows:

[TABLE]

These constraints are of the same form as the inequalities in Equation (8), for example the constraint $\widetilde{M}_{j}\geq 0$ can be written by choosing the normal vector $n_{i}$ to be the unit vector $-e_{j}^{nm}$ and $b_{i}$ to be the null vector $\mathbf{0}^{nm}$.

We note that there is a special case when the possible values of the parameters of the polynomial system are given as intervals, and the polyhedron $\mathcal{P}$ is a cuboid.

It is possible to define a set $L$ of finitely many additional linear constraints on the variables to characterize a special property of the realizations, for example a set of reactions to be excluded, or mass conservation on a given level, see e.g. [45]. These constraints can affect not only the entries of the coefficient matrix $M$ but the Kirchhoff matrix of the realizations as well. If the Kirchhoff matrix $A_{k}$ of the realization is represented by the point $\widetilde{A_{k}}\in\mathbb{R}^{m^{2}-m}$ storing the off-diagonal elements, and $r$ is the number of constraints in the set $L$, then the equations can be written in the form

[TABLE]

where $\alpha_{i}\in\mathbb{R}^{nm}$, $\beta_{i}\in\mathbb{R}^{m^{2}-m}$ and $d_{i}\in\mathbb{R}$ hold for all $i\in\{1,\ldots,r\}$. These constraints do not change the general properties of the model, and as it will be shown in Section 2.3, it can be modelled as a linear programming problem.

In the case of the uncertain model, we will examine realizations assuming a fixed set of complexes. Therefore, the known parameters are the polyhedron $\mathcal{P}$, the set $L$ of constraints and the matrix $Y$. Hence a constrained uncertain kinetic system is referred to as the triple $[\mathcal{P},L,Y]$, but we will call it an uncertain kinetic system for brevity.

Definition 2.3**.**

A reaction network $(Y,A_{k})$ is called a realization of the uncertain kinetic system $[\mathcal{P},L,Y]$ if there exists a coefficient matrix $M\in\mathbb{R}^{n\times m}$ so that the equation $M=Y\cdot A_{k}$ holds, the point $\widetilde{M}$ is in the polyhedron $\mathcal{P}$ and the entries of the matrices $M$ and $A_{k}$ fulfil the set $L$ of constraints. Since the matrix $Y$ is fixed but the coefficients of the polynomial system can vary, this realization is referred to as the matrix pair $(M,A_{k})$.

2.3 Computational model

Assuming a fixed set of complexes, a realization $(M,A_{k})$ of an uncertain kinetic system $[\mathcal{P},L,Y]$ can be computed using a linear optimization framework.

In the constraint satisfaction or optimization model, the variables are the entries of the matrix $M$ and the off-diagonal entries of the matrix $A_{k}$. The constraints regarding the realizations of the uncertain model can be written as follows:

[TABLE]

Equations (11) ensure that the parameters of the dynamics correspond to a point of the polyhedron $\mathcal{P}$. Dynamical equivalence is defined by Equation (12), while Equations (13) and (14) are required for the Kirchhoff property of matrix $A_{k}$ to be fulfilled. Moreover, the constraints in the set $L$ can be written in the form of Equation (10).

The objective function of the optimization model can be defined according to the desired properties of the realization, for example in order to examine if the reaction $C_{i}\rightarrow C_{j}$ can be present in the reaction network or not, the objective can be defined as $\max[A_{k}]_{ji}$.

We apply the representation of realizations of the uncertain model as points of the Euclidean space $\mathbb{R}^{m^{2}-m+nm}$. The coordinates with indices $i\in\{1,\ldots,m^{2}-m\}$ characterize the Kirchhoff matrix of the realization and the remaining coordinates $j\in\{m^{2}-m+1,\ldots m^{2}-m+nm\}$ define the coefficient matrix $M$ of the polynomial system. Due to the linearity of the constraints in the computational model, the set of possible realizations of an uncertain kinetic system $[\mathcal{P},L,Y]$ is a convex bounded polyhedron denoted by $\mathcal{Q}$.

3 Structural analysis of realizations of the uncertain

kinetic model

In this section we summarize some of the special structural properties of the realizations of an uncertain kinetic system $[\mathcal{P},L,Y]$.

3.1 Superstructure property of the dense realizations

A dynamically equivalent or linearly conjugate realization of a kinetic system with a fixed set of complexes having maximal or minimal number of reactions is called dense or sparse realization, respectively [33, 20]. It is known that for any kinetic system there might be several different sparse realizations, however, the dense realization is structurally unique and it defines a superstructure among all realizations, see [19].

The directed graph $G(V,E)$ is called a superstructure with respect to a set $\mathcal{G}$ of directed graphs with labelled vertices, if it contains every graph in the set $\mathcal{G}$ as subgraph, and it is minimal under inclusion. By the definition it follows that for any set $\mathcal{G}$ there exists a superstructure graph and it is unique.

In the case of dynamical equivalent and linearly conjugate realizations of kinetic systems the superstructure is the reaction graph of a dense realization, that contains all the reaction graphs representing realizations of the kinetic system as subgraphs, not considering the edge weights. This means that the set of reactions that take place in any of the realizations is the same as the set of reactions in the dense realization.

Dense and sparse realizations can be introduced in the case of the uncertain model as well, that are useful during the structural analysis.

Definition 3.1**.**

A realization $(M,A_{k})$ of the uncertain kinetic system $[\mathcal{P},L,Y]$ is called a dense (sparse) realization if it has maximal (minimal) number of reactions.

It can be proved that the superstructure property holds for uncertain kinetic systems as well, and the proof is based on the same idea as in the non-uncertain case, see [44].

Proposition 3.2**.**

A dense realization $(M,A_{k})$ of an uncertain kinetic system $[\mathcal{P},L,Y]$ determines a superstructure among all realizations of the model.

Proof.

If the point $D$ in the polyhedron $\mathcal{Q}$ of possible realizations represents a dense realization, then the superstructure property is equivalent to the property that any coordinate with index $i\in\{1,\ldots,m^{2}-m\}$ of an arbitrary point in $\mathcal{Q}$ can be positive only if the same coordinate of $D$ is positive. Let us assume by contradiction that there is another realization $R\in\mathcal{Q}$ so that there is an index $j\in\{1,\ldots,m^{2}-m\}$ for which $D_{j}=0$ and $R_{j}>0$ hold.

Since the polyhedron $\mathcal{Q}$ is closed under convex combination, the point

[TABLE]

is also in $\mathcal{Q}$. The coordinates with indices of the set $\{1,\ldots,m^{2}-m\}$ of all the points in $\mathcal{Q}$ are nonnegative, therefore such a coordinate of the convex combination is positive if the corresponding coordinate of $D$ or $R$ is positive. Consequently, $T$ has more positive coordinates with indices $j\in\{1,\ldots,m^{2}-m\}$ than the dense realization does, which is a contradiction. ∎

It follows from Proposition 3.2 that the structure of the dense realization is unique. If there were two different dense realizations, then the reaction graphs representing them would contain each other as subgraphs, which implies that these graphs are structurally identical.

The dense and sparse realizations are useful for checking the structural uniqueness of the uncertain model.

Proposition 3.3**.**

The dense and sparse realizations of an uncertain kinetic system $[\mathcal{P},L,Y]$ have the same number of reactions if and only if all realizations of the model are structurally identical.

Proof.

According to the definitions if in the dense and sparse realizations there is the same number of reactions, then in all realizations there must be the same number of reactions. Since the structure of the dense realization is unique, there cannot be two realizations with the maximal number of reactions but different structures, therefore all realizations must be structurally identical to the dense realization.

The converse statement is trivial: If all the realizations of the model are structurally identical, then the dense and sparse realizations must have identical structures, too. ∎

3.2 Polynomial-time algorithm to determine dense realizations

A dense realization of the uncertain kinetic system can be computed by the application of a recursive polynomial-time algorithm. The basic principle of the method is similar to the one presented in [44]: To each reaction a realization is assigned where the reaction takes place, if it is possible. In general, the same realization can be assigned to several reactions. Therefore, there is no need to perform a separate computation step for each reaction. The convex combination of the assigned realizations is also a realization of the uncertain model. If all the coefficients of the convex combination are positive then all reactions that take place in any of the assigned realizations are present in the convex combination as well. Consequently, the obtained realization represents a dense realization, where all reactions are present that are possible.

The computation can be performed in polynomial time since it requires at most $m^{2}-m$ steps of LP optimization and some minor computation.

Remark 3.4**.**

It follows from the operation of the algorithm that if there are at least two realizations assigned to reactions as defined, then there are infinitely many dense realizations, since at least one coefficient of the convex combination can be chosen arbitrarily from the interval $(0,1)$.

In the algorithm the assigned realizations are represented as points in $\mathbb{R}^{m^{2}-m+nm}$ and are determined using the following procedure:

FindPositive$([\mathcal{P},L,Y],H)$ returns a pair $(R,B)$. The point $R\in\mathcal{Q}$ represents a realization of the uncertain model $[\mathcal{P},L,Y]$ for which the value of the objective function $\sum_{j\in H}R_{j}$ considering a set $H\subseteq\{1,\ldots,m^{2}-m\}$ of indices is maximal. The other returned object is a set $B$ of indices where $k\in B$ if and only if $Q_{k}>0$. If there is no realization fulfilling the constraints then the pair $(\mathbf{0},\emptyset)$ is returned.

In the algorithm we apply the arithmetic mean as convex combination, i.e. if the number of the assigned realizations is $k$ then all the coefficients of the convex combination are $\frac{1}{k}$.

Proposition 3.5**.**

The realization returned by Algorithm 1 is a dense realization of the uncertain kinetic system.

Proof.

Since the set of all possible solutions can be represented as a convex polyhedron, the point $Result$ computed as the convex combination of realizations is indeed a realization of the uncertain kinetic system $[\mathcal{P},L,Y]$. Let us assume by contradiction that the returned point $Result$ does not represent the dense realization. Then there is a reaction $(C_{i},C_{j})$ which is present in the dense realization but it does not take place in $Result$. By the operation of the algorithm it follows that there must be a realization assigned to the reaction $(C_{i},C_{j})$, consequently this reaction takes place in the realization computed as the convex combination of the assigned realizations as well. This is a contradiction. ∎

3.3 Core reactions of uncertain models

A reaction is called core reaction of a kinetic system if it is present in every realization of the kinetic system [34]. It is possible that there are no core reactions, but there can be several of them as well. If all the realizations are structurally identical, then by Proposition 3.3 it follows that each reaction is a core reaction. The notion of core reactions can be extended to the case of uncertain models in a straightforward way.

Definition 3.6**.**

A reaction $C_{i}\rightarrow C_{j}$ is called a core reaction of the uncertain kinetic system $[\mathcal{P},L,Y]$ if it is present in each realization of the model, considering all possible coefficient matrices $M$ for which $\widetilde{M}\in\mathcal{P}$ holds.

Let $[\mathcal{P},L,Y]$ and $[\mathcal{P^{\prime}},L,Y]$ be two uncertain kinetic systems considering the same sets of complexes and additional linear constraints so that the polyhedron $\mathcal{P^{\prime}}$ is a subset of $\mathcal{P}$. If the sets of core reactions in the models are denoted as $C_{\mathcal{P}}$ and $C_{\mathcal{P^{\prime}}}$, respectively, then $C_{\mathcal{P}}\subseteq C_{\mathcal{P^{\prime}}}$ must hold. This property holds even if $\mathcal{P^{\prime}}$ is a single point in $\mathbb{R}^{nm}$ and $[\mathcal{P^{\prime}},L,Y]$ is a kinetic system defined as an uncertain kinetic system.

The set of core reactions of an uncertain kinetic system can be computed using a polynomial-time algorithm. This method has been first published in [39] for a special case, where the coefficients of the polynomial system have to be in predefined intervals, therefore the polyhedron $\mathcal{P}$ is a cuboid. Since the model applies only the property that all the constraints characterizing the model are linear, it can be applied without any modification to uncertain kinetic systems as well.

The question whether a certain reaction is a core reaction of a kinetic model or not, can be answered by solving a linear optimization problem. If this question has to be decided for all possible reactions, the computation can be done more effectively than doing separate optimization steps for every reaction. The idea is to minimize the sum of variables representing the off-diagonal entries of the Kirchhoff matrix. Generally, several variables in the minimized sum are zero in the computed realization, which means that the reactions corresponding to these variables are not core reactions. This step is repeated with the remaining set of variables until the computation does not return any non-core reactions. Finally, the remaining variables need to be checked one-by-one.

In the algorithm we refer to sets of indices corresponding to the off-diagonal entries of the Kirchhoff matrix $A_{k}$ by their characteristic vectors. The set $B\subseteq\{1,\ldots,m^{2}-m\}$ represented by the vector $b\in\{0,1\}^{m^{2}-m}$, which is defined as

[TABLE]

The procedure applied during the computation is more formally the following:

FindNonCore$([\mathcal{P},L,Y],b)$ computes a realization of the uncertain kinetic system $[\mathcal{P},L,Y]$ represented as a point $R\in\mathbb{R}^{m^{2}-m+nm}$, for which the sum of the coordinates with indices in the set $B\in\{1,\ldots,m^{2}-m\}$ is minimal. The procedure returns the vector $c$, the characteristic vector of the set $C$ which contains the indices corresponding to zero entries of the Kirchhoff matrix of the realization $R$, i.e. $C\subseteq\{1,\ldots,m^{2}-m\}$ and $\left[i\in C\ \Longleftrightarrow\ R_{i}=0\right]$.

We also need to utilize some operations on the sets represented by their characteristic vectors:

$b*c$ represents the set $B\cap C$, i.e. it is an element-wise ‘logical and’

$\overline{c}$ represents the complement of the set $C$, i.e. it is an element-wise negation.

Proposition 3.7**.**

Algorithm 2* computes the set of core reactions of an uncertain kinetic system $[\mathcal{P},L,Y]$ in polynomial time.*

Proof.

Let us assume by contradiction that the algorithm does not return the proper set of core reactions. There can be two different types of error:

$a)$ Let us assume that there is an index $i$ for which the corresponding reaction is a core reaction, but according to the algorithm it is not. In this case there must be a realization $R$ computed by the algorithm so that $R_{i}$ is zero. This is a contradiction.

$b)$ Let us assume that there is an index $j$ for which the corresponding reaction is not a core reaction but the algorithm returns the opposite answer. Consequently, after the while loop of the computation (from line 8) the coordinate $b_{j}$ must be equal to $1$. Then the remaining possible core reactions are examined one by one, therefore the procedure FindNonCore$([\mathcal{P},L,Y],e_{j}^{m^{2}-m})$ is also applied. According to the assumption the realization $R$ computed by the procedure must be so that $R_{j}$ is zero, which also yields a contradiction.

The computation according to the algorithm can be performed in polynomial time, since it requires the solution of at most $m^{2}-m$ LP optimization problems and some additional minor computation steps. ∎

4 Algorithm to determine all possible reaction graph structures of uncertain models

In this section we introduce an algorithm for computing all possible reaction graph structures of an uncertain kinetic system $[\mathcal{P},L,Y]$. The proposed method is an improved version of the algorithm published in [43], where all the optimization steps can be done parallelly. We also give a proof of the correctness of the presented method. Before presenting the pseudocode of the algorithm, we give a brief explanation of its data structures and operating principles.

We represent reaction graph structures by binary sequences, where each entry encodes the presence or lack of a reaction. During the algorithm, all data (i.e. the Kirchhoff and the coefficient matrices) of the realizations are computed, but only the binary sequences encoding the directed graph structures are stored and returned as results.

According to the superstructure property described in Proposition 3.2, only the reactions belonging to the dense realization need representation and storage. Moreover, if there are core reactions as well, then the coordinates corresponding to these can also be omitted. Both sets can be computed in polynomial time as it has been presented in Sections 3.2 and 3.3.

Let us refer to the set of reactions in the dense realization and the set of core reactions in the uncertain kinetic system $[\mathcal{P},L,Y]$ as $D_{\mathcal{P}}$ and $C_{\mathcal{P}}$, respectively. Then a realization of the uncertain model $[\mathcal{P},L,Y]$ can be represented by a binary sequence $R$ of length $z$, where $z$ is the size of the set $D_{\mathcal{P}}\setminus C_{\mathcal{P}}$ of non-core reactions in the dense realization. To define the binary sequence $R$ it is necessary to fix an ordering on the set of non-core reactions. The coordinate $R_{i}$ is equal to $1$ if and only if the $i$th non-core reaction is present in the realization, otherwise it is zero.

It is easy to see that knowing its structure, a realization can be determined in polynomial time: For each reaction $C_{i}\rightarrow C_{j}$ which is known not to be present in the realization the constraint $[A_{k}]_{ji}=0$ needs to be added to the constraint set $L$, and a dense realization of the (constrained) model has to be computed. Since it is known that there exists a realization where all non-excluded reactions take place, all of them have to be present in the computed constrained dense realization, consequently it will have exactly the prescribed structure.

During the computation the initial substrings of the binary sequences have a special role. Therefore, for all $k\in\{1,\ldots z\}$ a special equivalence relation $=_{k}$ is defined on the binary sequences. We say that $R=_{k}W$ holds if for all $i\in\{1,\ldots k\}$ the coordinate $R_{i}$ is equal to $W_{i}$. The equivalence class of the relation $=_{k}$ that contains the sequence $R$ as a representative is referred to as $C_{k}(R)$. (We note that in general there are several representatives of an equivalence class.) The elements of an equivalence class $C_{k}(R)$ can be characterized by a set of linear constraints added to the model. According to this property and Proposition 3.2, the dense realization in $C_{k}(R)$ determines a superstructure among all the realizations in the same set. The procedure FindRealization applied during the algorithm computes dense realizations of the uncertain model determined by the initial substrings. A realization is referred to as a pair $(R,k)$ if the corresponding realization represents the dense realization in $C_{k}(R)$. The realizations represented by such pairs get stored for some time in a stack $S$, the command ‘push $(R,k)$ into $S$’ puts the pair $(R,k)$ into the stack and ‘pop from $S$’ takes a pair out of the stack and returns it. The number of elements in the stack $S$ is denoted by $size(S)$.

The result of the entire computation is collected in a binary array called $Exist$, where all the computed graph structures are stored. The indices of the elements are the sequences as binary numbers, and the value of element $Exist[R]$ is equal to $1$ if and only if a realization with the structure encoded by $R$ has been found.

Considering the data structures, the main difference between the proposed method and the algorithm presented in [43] is that the sequences encoding the reaction graph structures are stored in only one stack in our current solution. Furthermore, the optimization steps using the sequences popped from this stack can be run in parallel. However, in this case the use of the binary array $Exist$ is necessary.

Within the algorithm we repeatedly apply two subroutines:

FindRealization$((R,k),i)$ computes a dense realization of the uncertain kinetic system $[\mathcal{P},L,Y]$, for which the representing binary sequence $W$ is in $C_{k}(R)$, and for every index $j\in\{k+1,\ldots,i\}$ the coordinate $W_{j}$ is zero. It is possible that among the first $k$ coordinates there are more zeros than required, therefore the computed sequence $W$ is compared to the sequence $R$. The procedure returns the sequence $W$ only if $W=_{k}R$ holds, otherwise $-1$ is returned. If the optimization task is infeasible then the returned object is also $-1$.

FindNextOne$((R,k))$ returns the smallest index $i$ for which $k<i$ and $R_{i}=1$ hold. If there is no such index, i.e. $R_{j}$ is zero for all $k<j$, then it returns $z+1$, where we recall that $z$ is the length of the sequences that encode the graph structures.

Let the sequence $D=\mathbf{1}$ represent the dense realization. Then the pseudocode of the algorithm for computing all possible graph structures can be given as follows.

Using the description of the algorithm, we can give formal results about its main properties.

Proposition 4.1**.**

Algorithm 3* computes all possible reaction graph structures representing realizations of an uncertain kinetic system $[\mathcal{P},L,Y]$.*

Proof.

Let us assume by contradiction that there is a realization of the uncertain kinetic system $[\mathcal{P},L,Y]$ represented by the sequence $V$ which is not returned by Algorithm 3. Let $R$ be another sequence that was stored in the stack $S$ as $(R,p)$ at some point during the computation, for which $V=_{p}R$ holds and $p$ is the greatest such number. If $p=0$ then $D$ is suitable to be $R$, and by the operation of the algorithm it follows that $p<z$ holds. (If $p$ were equal to $z$, then $V$ would be equivalent to $R$ which is a contradiction.)

There is a point during the computation when $(R,p)$ is popped out from the stack $S$. Let us assume that $\text{FindNextOne}(R,p)$ returns $i$ and $\text{FindNextOne}(V,p)$ returns $j$. In this case $i\leq j$ must hold since $R$ represents the superstructure in $C_{p}(R)$ and if $i$ were equal to $j$ then $p$ would not be maximal.

For the examination of sequence $R$, the procedure $\text{FindRealization}((R,p),i)$ is applied first (line 10), and it must return a valid sequence $W_{1}$, since its constrains are fulfilled by the realization $V$ as well. If $\text{FindNextOne}(W_{1},p)$ is $j_{1}$ then $j_{1}\leq j$ must hold, since $W_{1}$ represents the superstructure in $C_{i}(W_{1})$ and $V$ is also in $C_{i}(W_{1})$. If $j_{1}$ was equal to $j$ then $p$ would not be maximal. Otherwise, the computation can be continued by calling the procedure $\text{FindRealization}((R,p),j_{1})$. It must return a valid sequence $W_{2}$ for which we get that $\text{FindNextOne}(W_{2},p)=j_{2}\leq j$ holds by applying similar reasoning as above.

These steps must lead to contradiction either by $p$ not being maximal or by creating an infinite increasing sequence of integers that has an upper bound.

It follows that every possible reaction graph structure that represents a realization of the uncertain kinetic system $[\mathcal{P},L,Y]$ is returned by the algorithm. ∎

Remark 4.2**.**

Since the calculations of procedure FindRealization$((R,k),i)$ are independent of the results of previous calls of the same procedure, the order of the calls is irrelevant regarding the result of the entire computation.

Remark 4.3**.**

The proof of Proposition 3.2 in [43] can be applied for verifying the property that during the computation according to Algorithm 3 every reaction graph structure is returned only once.

Remark 4.4**.**

We can also give an upper bound to the number of required optimization steps by considering the realizations $(R,k)$ regarding $k$. For all $k$ the number of possible realizations $R$ stored in the stack $S$ is at most $2^{k}$. When such a realization is popped from the stack the required optimization steps is at most $z-k$. Consequently, a rough upper bound to the number of optimization steps required during Algorithm 3 can be given as $\sum\limits_{k=0}^{z-1}2^{k}(z-k)$.

5 Illustrative examples

In this section we demonstrate the operation of the algorithms presented in this paper on two examples in case of different degrees and types of uncertainties, and even in the case of additional linear constraints.

5.1 Example 1: a simple kinetic system

The model that serves as a basis for this example was presented previously in [36, 45]. The uncertain model is generated using the kinetic system

[TABLE]

where $c_{1},c_{2}>0$. We consider realizations on a fixed set $\mathcal{C}=\{C_{1},C_{2},C_{3}\}$ of complexes, where the complexes $C_{1}=3X_{2}$, $C_{2}=3X_{1}$, $C_{3}=2X_{1}+X_{2}$ are formed of the species $X_{1}$ and $X_{2}$. It follows that the characterizing matrices $Y$ and $M$ of the kinetic system referred to as $[M,Y]$ are

$Y=\begin{bmatrix}0&3&2\\ 3&0&1\end{bmatrix}\qquad M=\begin{bmatrix}[r]3c_{1}&-c_{2}&0\\ -3c_{1}&c_{2}&0\end{bmatrix}$

During the numerical computations the parameter values $c_{1}=1$ and $c_{2}=2$ were applied.

*A. Uncertainty defined by independent intervals

*This model represents a special case in the class of uncertain kinetic systems defined in Section 2.2, since the possible values of every coefficient of the kinetic system are determined by independent upper and lower bounds that are defined as relative distances. Let us represent the entry $[M]_{ij}$ of the coefficient matrix $M$ by the coordinate $\widetilde{M}_{l}$ of the point $\widetilde{M}\in\mathbb{R}^{6}$. Moreover, let the relative distances of the upper and lower bounds of $\widetilde{M}_{l}$ be given by the real constants $\gamma_{l}$ and $\rho_{l}$ from the interval $[0,1]$, respectively. Then the equations defining the polyhedron $\mathcal{P}_{A}\subset\mathbb{R}^{6}$ of the uncertain parameters can be written in terms of the coordinates $\widetilde{M}_{l}$ as

[TABLE]

In the examined uncertain kinetic system $[\mathcal{P}_{A},L,Y]$ no additional linear constraints are considered, i.e. $L=\emptyset$.

In [45] all possible reaction graphs – with the indication of the reaction rate constants defined as functions of the parameters $c_{1}$ and $c_{2}$ – representing dynamically equivalent realizations of the kinetic system $[M,Y]$ have been presented. Obviously, these structures must appear among the realizations of the uncertain kinetic model $[\mathcal{P}_{A},\emptyset,Y]$ as well, but there might be additional possible structures among the realizations of the uncertain kinetic system.

Interestingly, the result of the computation was that in the case of any degree of uncertainty ($\gamma_{l},\rho_{l}\in[0,1)$ for all $l\in\{1,\ldots 6\}$), the sets of possible reaction graph structures of the uncertain model $[\mathcal{P}_{A},\emptyset,Y]$ and that of the non-uncertain system $[M,Y]$ are identical. This result might be contrary to expectations, but for this small example it is easy to prove that the obtained graph structures are indeed correct for all positive values of the parameters $c_{1}$ and $c_{2}$. For this, we divide the computation into smaller steps.

It has been shown in [43] that in the case of dynamically equivalent realizations the computation can be done column-wise (since the $j$th column of matrix $A_{k}$ depends only on the $j$th column of matrix $M$). These computations can be performed separately, and all the possible reaction graph structures can be constructed by choosing a column structure for every index $j\in\{1,\ldots,m\}$ and building the Kirchhoff matrix $A_{k}$ of the realization from them. Consequently, if in the case of the $j$th column the number of different structures is $p_{j}$, then the number of structurally different realizations is $\prod\limits_{j=1}^{m}p_{j}$.

First the original kinetic system $[M,Y]$ is examined. To make the notations less complicated, the entries of the Kirchhoff matrix are denoted by the corresponding reaction rate coefficients, i.e. $[A_{k}]_{ij}=k_{ji}$ for all $i,j\in\{1,2,3\}$, $i\neq j$.

In the case of the first column we get:

[TABLE]

It can be seen that for every positive value of the parameter $c_{1}$ the two corresponding reaction rates can realize 3 of the $2^{2}=4$ possible structurally different solutions. Both can be positive, or either one can be positive while the other one is zero. (Possible outcomes are for example: $k_{12}=\frac{1}{2}c_{1},k_{13}=\frac{3}{4}c_{1}$ or $k_{12}=0,k_{13}=\frac{3}{2}c_{1}$ or $k_{12}=c_{1},k_{13}=0$.) The fourth case, when both $k_{12}$ and $k_{13}$ are zero is possible only when $[M]_{.1}=[0\leavevmode\nobreak\ \leavevmode\nobreak\ 0]^{\top}$, which requires the corresponding parameters of uncertainty $\rho_{i}$ to be at least one.

In the case of the second column, 3 of the 4 possible outcomes can be realized and a similar reasoning can be applied:

[TABLE]

In the third column there is no uncertainty because there are only zero entries in $[M]_{.3}$. Consequently, in the case of $[A_{k}]_{.3}$ only 2 solutions are possible. The two corresponding reactions can either be both present or both missing.

[TABLE]

It follows from the above computations that the number of possible reaction graph structures is $3\cdot 3\cdot 2=18$, and the generated structures are identical to the ones presented in [45]. This number could be larger only if all the reaction rates in the first or second column of $A_{k}$ can be zero, but this requires the entries in the corresponding column $[M]_{.1}$ or $[M]_{.2}$ to be zero.

*B. Uncertainty defined as a general polyhedron

*Now we examine the uncertain kinetic system that was also generated from the kinetic system $[M,Y]$, but the set $\mathcal{P}_{B}$ of possible coefficients is defined as a more general polyhedron.

If the matrix $M$ of coefficients is represented by the point $\widetilde{M}\in\mathbb{R}^{6}$ so that $\widetilde{M_{\ }}^{\top}=[M_{11},M_{12},M_{13},M_{21},M_{22},M_{23}]$, then let the equations determining the polyhedron $\mathcal{P}_{B}$ be the following:

[TABLE]

In this case, again, no additional linear constraints are considered in the uncertain model, i.e. we examine the uncertain model $[\mathcal{P}_{B},\emptyset,Y]$. The computation of all possible reaction graph structures shows that in addition to the structures realizing the non-uncertain kinetic system $[M,Y]$, there are 6 more possible structures, presented in Figure 1.

It can be seen that the point $\widetilde{M}_{1}^{\top}=[3,-2,0,-3,2,0]$ corresponding to the original kinetic system is in the polyhedron $\mathcal{P}_{B}$, therefore the 18 structures determined by its realizations must be among the realizations of the uncertain kinetic system. Then we can apply a reasoning similar to that in Section 5.1.$A$. Since the entries in column $[M]_{.3}$ are all zero in every point of $\mathcal{P}_{B}$, only the two outcomes that appear in the case of the original kinetic system $[M,Y]$ are possible in the case of this column. The uncertain model can have more realizations than the original kinetic system only if all the reaction rates in at least one of the columns $[A_{k}]_{.1}$ or $[A_{k}]_{.2}$ can be zero. This is possible only if all the entries in $[M]_{.1}$ or $[M]_{.2}$ are zero. From the constraints of the polyhedron $\mathcal{P}_{B}$ it follows that $[M]_{11}\geq 1$, consequently the column $[M]_{.1}$ cannot be zero. But $[M]_{.2}$ can have only zero entries, for example the point $\widetilde{M}_{2}=[3,0,0,-3,0,0]^{\top}\in\mathcal{P}_{B}$ satisfies this property. For the columns of the matrices $M$ and $M_{2}$ the following hold: $[M_{2}]_{.1}=[M]_{.1}$ and $[M_{2}]_{.3}=[M]_{.3}$. Therefore, for the first and third columns of $A_{k}$ there are 3 and 2 possible outcomes, respectively. Since in the case of the second column there is one additional possible outcome, the number of further reaction graph structures (compared to the original kinetic system $[M,Y]$) is $3\cdot 2=6$. It is easy to see that these are exactly the ones presented in Figure 1 with the indicated reaction rate coefficients for an arbitrary $p>0$.

5.2 Example 2: G-protein network

The G-protein (guanine nucleotide-binding protein) cycle has a key role in several intracellular signalling transduction pathways. The G-protein located on the intracellular surface of the cell membrane is activated by the binding of specific ligand molecules to the G-protein coupled receptor of the extracellular membranes surface. The activated G-protein dissociates to different subunits which take part in intracellular signalling pathways. After the termination of the signalling mechanisms, the subunits become inactive and bind into each other [27].

We examined the structural properties of the yeast G-protein cycle using the model published in [42]. The model involves a so-called heterotrimeric G-protein containing three different subunits. In response to the extracellular ligand binding, the protein dissociates to G-$\alpha$ and G-$\beta\gamma$ subunits, where the active and inactive forms of the G-$\alpha$ subunit can also be distinguished.

The reaction network model involves the following species: $R$ and $L$ represent the receptor and the corresponding ligand, respectively, $RL$ refers to the ligand-bound receptor, $G$ is the G-protein located on the intracellular membrane surface, $G_{a}$ and $G_{d}$ denote the active and the inactive forms of the G-$\alpha$ subunit and $G_{bg}$ is the G-$\beta\gamma$ subunit.

The model can be characterized as a chemical reaction network $(Y,A_{k})$, where the structures of the complexes and the reactions are defined by the complex composition matrix $Y\in\mathbb{R}^{7\times 10}$ and the Kirchhoff matrix $A_{k}\in\mathbb{R}^{10\times 10}$ as follows:

[TABLE]

The kinetic system that is realized by the model is $\dot{x}=M\cdot\psi^{Y}=Y\cdot A_{k}\cdot\psi^{Y}$, i.e. $M=Y\cdot A_{k}\in\mathbb{R}^{7\times 10}$. The reaction graph structure of the G-protein model can be seen in Figure 2 with the indication of the linkage classes. (The linkage classes are the undirected connected components of the reaction graph.)

The computation of all possible reaction graph structures and the solution of the linear equations shows that the heterotrimeric G-protein cycle with the given parametrization is not just structurally but also parametrically unique. Thus the prescribed dynamics without uncertainty cannot be realized by any other set of reactions or different reaction rate coefficients using the given set of complexes.

A. Uncertainty defined with independent relative distance intervals

We have examined the uncertain kinetic systems defined by relative parameter uncertainty as it was presented in Section 5.1.

First we examined the uncertain model $[\mathcal{P}_{0.1},\emptyset,Y]$, where the uncertainty coefficients $\gamma_{l}$ and $\rho_{l}$ for all $l\in\{1,\ldots,70\}$ are $0.1$ and there are no additional linear constraints in the model. By computing all possible reaction graph structures and the set of core reactions of this uncertain kinetic system, we obtained that all the reactions in the original G-protein cycle are core reactions. Moreover, in the dense realization there are 10 further reactions, and these can be present in the realization independently of each other. Consequently, the total number of different graph structures is $2^{10}=1024$. Figure 3 shows the number of possible reaction graph structures with different number of reactions. The dense realization for this case is shown in Figure 4.

If we increase the relative uncertainty to 0.2, we obtain the uncertain kinetic system $[\mathcal{P}_{0.2},\emptyset,Y]$ with $\gamma_{l}=\rho_{l}=0.2$ for $i=1,\dots,70$. In this case, the reaction $RL\ \rightarrow 0$ is no longer a core reaction, and it can also be added or removed independently of all other reactions (which remain independent of each other). Therefore, the number of possible structures becomes $2^{11}=2048$.

B. Constrained uncertain model

We have also examined the possible structures in the case of constrained uncertain models. The set $L_{1}$ of constraints prohibits every reaction between different linkage classes. It can be seen in Figure 5 that the dense realization of the uncertain kinetic system $[\mathcal{P}_{0.1},L_{1},Y]$ has 3 reactions that are exactly the ones that are present in the dense realization of $[\mathcal{P}_{0.1},\emptyset,Y]$ and do not connect different linkage classes. These reactions are independent of each other, therefore the number of structurally different realizations is $2^{3}=8$ in the case of the uncertain kinetic system $[\mathcal{P}_{0.1},L_{1},Y]$ and $2^{4}=16$ for the model $[\mathcal{P}_{0.2},L_{1},Y]$. The sets of core reactions are the same as in the case of the unconstrained model for both degrees of uncertainty.

The independence of non-core reactions is a special property of the studied uncertain model. As a consequence of this and the superstructure property of the dense realization, the dense realization of the constrained model will contain each reaction of the unconstrained model that is not excluded by the constraints.

We emphasize that the dense realizations in the above example contain all mathematically possible reactions that can be compatible with the studied uncertain models. If, using prior knowledge, the biologically non-plausible reactions are excluded and/or certain relations between model parameters are ensured via linear constraints, then the described methodology is still suitable to check the structural uniqueness of the resulting uncertain kinetic model.

6 Conclusion

The set of reaction graph structures realizing uncertain kinetic models was studied in this paper. For this, an uncertain polynomial model class was introduced, where the coefficients of monomials belong to a polytopic set. Thus, an uncertain kinetic model includes a set of kinetic ordinary differential equations. Using the convexity of the parameter set, it was proved that the unweighted dense reaction graph containing the maximum number of reactions corresponding to an uncertain model, forms a superstructure among the possible realizations assuming a fixed complex set. This means that any unweighted reaction graph realizing any kinetic ODE within an uncertain model is a subgraph of the unweighted directed graph of the dense realization.

To search through the possible graph structures, an optimization-based computational model was introduced, where the decision variables are the reaction rate coefficients, and the entries of the monomial coefficient matrix. It was shown that the dense realization can be computed in polynomial time using linear programming steps. An algorithm was proposed to compute those ‘invariant’ reactions (called core reactions) of uncertain models, that are present in any realization of the uncertain model. Most importantly, an algorithm with correctness proof was also proposed in the paper for enumerating all possible reaction graph structures for an uncertain kinetic model.

The theoretical results and proposed algorithms were illustrated on two examples. The examples show that the proposed approach is suitable for the structural uniqueness analysis of uncertain kinetic models.

Acknowledgements

This project was developed with the support of the PhD program of the Roska Tamás Doctoral School of Sciences and Technology, Faculty of Information Technology and Bionics, Pázmány Péter Catholic University, Budapest. The authors gratefully acknowledge the support of grants PPKE KAP-1.1-16-ITK, and K115694 of the National Research, Development and Innovation Office - NKFIH.