A computational approach to the structural analysis of uncertain kinetic systems

# A computational approach to the structural analysis of uncertain kinetic systems

## 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.

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

the set of real numbers
the set of nonnegative real numbers
the set of natural numbers
the set of matrices having entries from a set with rows and columns
the entry of matrix with row index and column index
the th column of matrix
the th coordinate of vector
the null vector in
a vector in with all coordinates equal to 1
a vector in for which the th coordinate is 1 and all the others are zero

### 2.1 Kinetic polynomial systems and their models

Nonnegative polynomial systems are defined in the following general form:

 ˙x=M⋅φ(x) (1)

where is a nonnegative valued function, is a coefficient matrix and 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 depending on the exponents of , 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: complexes: reactions:
For all , the reaction is represented by the ordered pair , and it is described by a nonnegative real number called reaction rate coefficient. The reaction is present in the reaction network if and only if is strictly positive.

The relation between species and complexes is described by the complex composition matrix , the columns of which correspond to the complexes, i.e.

 [Y]ij=αjii∈{1,…,n}, j∈{1,…,m} (2)

The presence of the reactions in the CRN is defined through the rate coefficients as the off-diagonal entries of the Kirchhoff matrix which is a Metzler compartmental matrix with zero column-sums. Its entries are defined as:

 [Ak]ij=⎧⎪⎨⎪⎩kji if i≠j−m∑l=1,l≠ikil if i=ji,j∈{1,…,m} (3)

According to this notation, the reaction takes place in the reaction network if and only if is positive, and implies that . Since a chemical reaction network is uniquely characterized by the matrices and , we refer to a CRN by the corresponding pair .

If mass action kinetics is assumed, the equations governing the dynamics of the concentrations of the species in the CRN defined by the function can be written in the form:

 ˙x=Y⋅Ak⋅ψY(x) (4)

where is the monomial function of the CRN with coordinate functions

 ψYj(x)=n∏i=1x[Y]iji,j∈{1,…,m} (5)

The nonnegative polynomial system (1) is called a kinetic system if there exists a reaction network so that its dynamics satisfies the equation [46]:

 M⋅φ(x)=Y⋅Ak⋅ψY(x) (6)

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 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 is equal to (and 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:

 M=Y⋅Ak (7)

Reaction networks have another representation, which is more suitable for illustrating the structural properties. It is a weighted directed graph 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 and correspond to the complexes and , respectively. Then there is a directed edge with weight if and only if the reaction takes place in the CRN.

### 2.2 Uncertain kinetic systems

For the uncertainty modelling, we assume that the monomial coefficients in matrix are constant but uncertain, and they belong to an 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 as a point denoted by in the Euclidean space . In the uncertain model it is assumed that the possible points are all the points of a closed convex polyhedron , which is defined as the intersection of halfspaces. The boundaries of the halfspaces are hyperplanes with normal vectors and constants . Applying these notations, the polyhedron can be described by a linear inequality system as

 P={˜M∈Rnm | ˜M⊤⋅ni≤bi, 1≤i≤q}} (8)

For the characterization of the polyhedron 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 as follows:

 [Y]ij=0 ⟹ [M]ij≥0,i∈{1,…,n}, j∈{1,…,m} (9)

These constraints are of the same form as the inequalities in Equation (8), for example the constraint can be written by choosing the normal vector to be the unit vector and to be the null vector .

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 is a cuboid.

It is possible to define a set 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 but the Kirchhoff matrix of the realizations as well. If the Kirchhoff matrix of the realization is represented by the point storing the off-diagonal elements, and is the number of constraints in the set , then the equations can be written in the form

 ˜M⊤⋅αi+˜Ak⊤⋅βi≤di (10)

where , and hold for all . 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 , the set of constraints and the matrix . Hence a constrained uncertain kinetic system is referred to as the triple , but we will call it an uncertain kinetic system for brevity.

###### Definition 2.3.

A reaction network is called a realization of the uncertain kinetic system if there exists a coefficient matrix so that the equation holds, the point is in the polyhedron and the entries of the matrices and fulfil the set of constraints. Since the matrix is fixed but the coefficients of the polynomial system can vary, this realization is referred to as the matrix pair .

### 2.3 Computational model

Assuming a fixed set of complexes, a realization of an uncertain kinetic system can be computed using a linear optimization framework.

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

 ˜M⊤⋅ni≤bi, i∈{1,…,p} (11) M=Y⋅Ak, (12) [Ak]ij≥0, i≠j, i,j∈{1,…,m} (13) m∑j=1[Ak]ij=0, j∈{1,…,m} (14)

Equations (11) ensure that the parameters of the dynamics correspond to a point of the polyhedron . Dynamical equivalence is defined by Equation (12), while Equations (13) and (14) are required for the Kirchhoff property of matrix to be fulfilled. Moreover, the constraints in the set 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 can be present in the reaction network or not, the objective can be defined as .

We apply the representation of realizations of the uncertain model as points of the Euclidean space . The coordinates with indices characterize the Kirchhoff matrix of the realization and the remaining coordinates define the coefficient matrix 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 is a convex bounded polyhedron denoted by .

## 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 .

### 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 is called a superstructure with respect to a set of directed graphs with labelled vertices, if it contains every graph in the set as subgraph, and it is minimal under inclusion. By the definition it follows that for any set 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 of the uncertain kinetic system 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 of an uncertain kinetic system determines a superstructure among all realizations of the model.

###### Proof.

If the point in the polyhedron of possible realizations represents a dense realization, then the superstructure property is equivalent to the property that any coordinate with index of an arbitrary point in can be positive only if the same coordinate of is positive. Let us assume by contradiction that there is another realization so that there is an index for which and hold.

Since the polyhedron is closed under convex combination, the point

 T=c⋅D+(1−c)⋅Rc∈(0,1)

is also in . The coordinates with indices of the set of all the points in are nonnegative, therefore such a coordinate of the convex combination is positive if the corresponding coordinate of or is positive. Consequently, has more positive coordinates with indices 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 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 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 .

In the algorithm the assigned realizations are represented as points in and are determined using the following procedure:

FindPositive returns a pair . The point represents a realization of the uncertain model for which the value of the objective function considering a set of indices is maximal. The other returned object is a set of indices where if and only if . If there is no realization fulfilling the constraints then the pair is returned.

In the algorithm we apply the arithmetic mean as convex combination, i.e. if the number of the assigned realizations is then all the coefficients of the convex combination are .

###### 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 computed as the convex combination of realizations is indeed a realization of the uncertain kinetic system . Let us assume by contradiction that the returned point does not represent the dense realization. Then there is a reaction which is present in the dense realization but it does not take place in . By the operation of the algorithm it follows that there must be a realization assigned to the reaction , 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 is called a core reaction of the uncertain kinetic system if it is present in each realization of the model, considering all possible coefficient matrices for which holds.

Let and be two uncertain kinetic systems considering the same sets of complexes and additional linear constraints so that the polyhedron is a subset of . If the sets of core reactions in the models are denoted as and , respectively, then must hold. This property holds even if is a single point in and 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 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 by their characteristic vectors. The set represented by the vector , which is defined as

 bi={1 if i∈B0 if i∉B (15)

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

FindNonCore computes a realization of the uncertain kinetic system represented as a point , for which the sum of the coordinates with indices in the set is minimal. The procedure returns the vector , the characteristic vector of the set which contains the indices corresponding to zero entries of the Kirchhoff matrix of the realization , i.e. and .

We also need to utilize some operations on the sets represented by their characteristic vectors: represents the set , i.e. it is an element-wise ‘logical and’ represents the complement of the set , i.e. it is an element-wise negation.

###### Proposition 3.7.

Algorithm 2 computes the set of core reactions of an uncertain kinetic system 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:

Let us assume that there is an index 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 computed by the algorithm so that is zero. This is a contradiction.

Let us assume that there is an index 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 must be equal to . Then the remaining possible core reactions are examined one by one, therefore the procedure FindNonCore is also applied. According to the assumption the realization computed by the procedure must be so that 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 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 . 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 as and , respectively. Then a realization of the uncertain model can be represented by a binary sequence of length , where is the size of the set of non-core reactions in the dense realization. To define the binary sequence it is necessary to fix an ordering on the set of non-core reactions. The coordinate is equal to if and only if the 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 which is known not to be present in the realization the constraint needs to be added to the constraint set , 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 a special equivalence relation is defined on the binary sequences. We say that holds if for all the coordinate is equal to . The equivalence class of the relation that contains the sequence as a representative is referred to as . (We note that in general there are several representatives of an equivalence class.) The elements of an equivalence class 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 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 if the corresponding realization represents the dense realization in . The realizations represented by such pairs get stored for some time in a stack , the command ‘push into ’ puts the pair into the stack and ‘pop from ’ takes a pair out of the stack and returns it. The number of elements in the stack is denoted by .

The result of the entire computation is collected in a binary array called , where all the computed graph structures are stored. The indices of the elements are the sequences as binary numbers, and the value of element is equal to if and only if a realization with the structure encoded by 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 is necessary.

Within the algorithm we repeatedly apply two subroutines:

FindRealization computes a dense realization of the uncertain kinetic system , for which the representing binary sequence is in , and for every index the coordinate is zero. It is possible that among the first coordinates there are more zeros than required, therefore the computed sequence is compared to the sequence . The procedure returns the sequence only if holds, otherwise is returned. If the optimization task is infeasible then the returned object is also .

FindNextOne returns the smallest index for which and hold. If there is no such index, i.e. is zero for all , then it returns , where we recall that is the length of the sequences that encode the graph structures.

Let the sequence 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 .

###### Proof.

Let us assume by contradiction that there is a realization of the uncertain kinetic system represented by the sequence which is not returned by Algorithm 3. Let be another sequence that was stored in the stack as at some point during the computation, for which holds and is the greatest such number. If then is suitable to be , and by the operation of the algorithm it follows that holds. (If were equal to , then would be equivalent to which is a contradiction.)

There is a point during the computation when is popped out from the stack . Let us assume that returns and returns . In this case must hold since represents the superstructure in and if were equal to then would not be maximal.

For the examination of sequence , the procedure is applied first (line 10), and it must return a valid sequence , since its constrains are fulfilled by the realization as well. If is then must hold, since represents the superstructure in and is also in . If was equal to then would not be maximal. Otherwise, the computation can be continued by calling the procedure . It must return a valid sequence for which we get that holds by applying similar reasoning as above.

These steps must lead to contradiction either by 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 is returned by the algorithm.

###### Remark 4.2.

Since the calculations of procedure FindRealization 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 regarding . For all the number of possible realizations stored in the stack is at most . When such a realization is popped from the stack the required optimization steps is at most . Consequently, a rough upper bound to the number of optimization steps required during Algorithm 3 can be given as .

## 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

 ˙x1=3c1⋅x32−c2⋅x31 ˙x2=−3c1⋅x32+c2⋅x31, (16)

where . We consider realizations on a fixed set of complexes, where the complexes , , are formed of the species and . It follows that the characterizing matrices and of the kinetic system referred to as are

During the numerical computations the parameter values and 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 of the coefficient matrix by the coordinate of the point . Moreover, let the relative distances of the upper and lower bounds of be given by the real constants and from the interval , respectively. Then the equations defining the polyhedron of the uncertain parameters can be written in terms of the coordinates as

 ˜M⊤⋅e6l≤(1+γl)⋅[M]ij (17) ˜M⊤⋅(−e6l)≤(ρl−1)⋅[M]ij (18)

In the examined uncertain kinetic system no additional linear constraints are considered, i.e. .

In [45] all possible reaction graphs – with the indication of the reaction rate constants defined as functions of the parameters and – representing dynamically equivalent realizations of the kinetic system have been presented. Obviously, these structures must appear among the realizations of the uncertain kinetic model 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 ( for all ), the sets of possible reaction graph structures of the uncertain model and that of the non-uncertain system 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 and . 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 th column of matrix depends only on the th column of matrix ). These computations can be performed separately, and all the possible reaction graph structures can be constructed by choosing a column structure for every index and building the Kirchhoff matrix of the realization from them. Consequently, if in the case of the th column the number of different structures is , then the number of structurally different realizations is .

First the original kinetic system is examined. To make the notations less complicated, the entries of the Kirchhoff matrix are denoted by the corresponding reaction rate coefficients, i.e. for all , .

In the case of the first column we get:

 Y⋅⎡⎢⎣−k12−k13k12k13⎤⎥⎦=[3c1−3c1]k12,k13∈R+⟹k12∈[0,c1], k13=32c1−32k12 (19)

It can be seen that for every positive value of the parameter the two corresponding reaction rates can realize 3 of the 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: or or .) The fourth case, when both and are zero is possible only when , which requires the corresponding parameters of uncertainty 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:

 Y⋅⎡⎢⎣k21−k21−k23k23⎤⎥⎦=[−c2c2]k21,k23∈R+⟹k21∈(0,c23),k23=c2−3k21 (20)

In the third column there is no uncertainty because there are only zero entries in . Consequently, in the case of only 2 solutions are possible. The two corresponding reactions can either be both present or both missing.

 Y⋅⎡⎢⎣k31k32−k31−k32⎤⎥⎦=[00]k31,k32∈R+⟹k31∈R+,k32=2k31 (21)

It follows from the above computations that the number of possible reaction graph structures is , 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 can be zero, but this requires the entries in the corresponding column or 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 , but the set of possible coefficients is defined as a more general polyhedron.

If the matrix of coefficients is represented by the point so that , then let the equations determining the polyhedron be the following:

 ˜M⊤⋅(−e61) ≤0 ˜M⊤⋅(−e65) ≤0 ˜M⊤⋅e63 =0 ˜M⊤⋅e66 =0 (22) ˜M⊤⋅[1,1,0,1,1,0]⊤ =0 ˜M⊤⋅[0,−1,0,−1,0,0]⊤ ≤7 ˜M⊤⋅[−1,0,0,0,1,0]⊤ ≤−1

In this case, again, no additional linear constraints are considered in the uncertain model, i.e. we examine the uncertain model . The computation of all possible reaction graph structures shows that in addition to the structures realizing the non-uncertain kinetic system , there are 6 more possible structures, presented in Figure 1.

It can be seen that the point corresponding to the original kinetic system is in the polyhedron , 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.. Since the entries in column are all zero in every point of , only the two outcomes that appear in the case of the original kinetic system 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 or can be zero. This is possible only if all the entries in or are zero. From the constraints of the polyhedron it follows that , consequently the column cannot be zero. But can have only zero entries, for example the point satisfies this property. For the columns of the matrices and the following hold: and . Therefore, for the first and third columns of 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 ) is . It is easy to see that these are exactly the ones presented in Figure 1 with the indicated reaction rate coefficients for an arbitrary .

### 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- and G- subunits, where the active and inactive forms of the G- subunit can also be distinguished.

The reaction network model involves the following species: and represent the receptor and the corresponding ligand, respectively, refers to the ligand-bound receptor, is the G-protein located on the intracellular membrane surface, and denote the active and the inactive forms of the G- subunit and is the G- subunit.

The model can be characterized as a chemical reaction network , where the structures of the complexes and the reactions are defined by the complex composition matrix and the Kirchhoff matrix as follows:

 Y =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣1010000000001000000001000010000001001000000010010000000001100000010010⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ Ak =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣−0.40000000040000−140.3220000000010−0.322000000000000000100000000−110000000000001100000000000000−0.010000000000.0100000000000−100000.440000000−4000⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

The kinetic system that is realized by the model is , i.e. . 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 , where the uncertainty coefficients and for all are 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 . 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 with for . In this case, the reaction 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 .

B. Constrained uncertain model

We have also examined the possible structures in the case of constrained uncertain models. The set 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 has 3 reactions that are exactly the ones that are present in the dense realization of and do not connect different linkage classes. These reactions are independent of each other, therefore the number of structurally different realizations is in the case of the uncertain kinetic system and for the model . 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.

### Footnotes

1. corresponding author

### References

1. D. F. Anderson. A proof of the Global Attractor Conjecture in the single linkage class case. SIAM Journal on Applied Mathematics, 71:1487–1508, 2011. http://arxiv.org/abs/1101.0761,.
2. V. Badri, M. J. Yazdanpanah, and M. S. Tavazoei. On stability and trajectory boundedness of Lotka–Volterra systems with polytopic uncertainty. IEEE Transactions on Automatic Control, pages available online, DOI 10.1109/TAC.2017.2663839, 2017.
3. S. Boyd, L. El-Ghaoui, E. Feron, and V. Balakrishnan. Linear Matrix Inequalities in Systems and Control Theory. SIAM Books, Philadelphia, PA, 1994.
4. W. Briggs. Uncertainty: The Soul of Modeling, Probability & Statistics. Springer, 2016.
5. V. Chellaboina, S. P. Bhat, W. M. Haddad, and D. S. Bernstein. Modeling and analysis of mass-action kinetics – nonnegativity, realizability, reducibility, and semistability. IEEE Control Systems Magazine, 29:60–78, 2009.
6. W. W. Chen, M. Niepel, and P. K. Sorger. Classic and contemporary approaches to modeling biochemical reactions. Geners & Development, 24:1861–1875, 2012.
7. O. T. Chis, J. R. Banga, and E. Balsa-Canto. Structural identifiability of systems biology models: a critical comparison of methods. PLOS One, 6(11):e27755, 2011.
8. O. T. Chis, A. F. Villaverde, J. R. Banga, and E. Balsa-Canto. On the relationship between sloppiness and identifiability. Mathematical Biosciences, 282:147–161, 2016.
9. G. Craciun. Toric differential inclusions and a proof of the global attractor conjecture. arXiv:1501.02860 [math.DS], January 2015.
10. G. Craciun and C. Pantea. Identifiability of chemical reaction networks. Journal of Mathematical Chemistry, 44:244–259, 2008.
11. M. Feinberg. Lectures on chemical reaction networks. Notes of lectures given at the Mathematics Research Center, University of Wisconsin, 1979.
12. M. Feinberg. Chemical reaction network structure and the stability of complex isothermal reactors - I. The deficiency zero and deficiency one theorems. Chemical Engineering Science, 42 (10):2229–2268, 1987.
13. T. V. Guy, M. Karny, and D. H. Wolpert, editors. Decision Making: Uncertainty, Imperfection, Deliberation and Scalability. Springer, 2015.
14. W. M. Haddad, VS. Chellaboina, and Q. Hui. Nonnegative and Compartmental Dynamical Systems. Princeton University Press, 2010.
15. G. W. Harrison. compartmental models with uncertain flow rates. Mathematical Biosciences, 43:131–139, 1979.
16. F. Horn and R. Jackson. General mass action kinetics. Archive for Rational Mechanics and Analysis, 47:81–116, 1972.
17. V. Hárs and J. Tóth. On the inverse problem of reaction kinetics. In M. Farkas and L. Hatvani, editors, Qualitative Theory of Differential Equations, volume 30 of Coll. Math. Soc. J. Bolyai, pages 363–379. North-Holland, Amsterdam, 1981.
18. M. D. Johnston and D. Siegel. Linear conjugacy of chemical reaction networks. Journal of Mathematical Chemistry, 49:1263–1282, 2011.
19. M. D. Johnston, D. Siegel, and G. Szederkényi. Dynamical equivalence and linear conjugacy of chemical reaction networks: new results and methods. MATCH Commun. Math. Comput. Chem., 68:443–468, 2012.
20. M. D. Johnston, D. Siegel, and G. Szederkényi. A linear programming approach to weak reversibility and linear conjugacy of chemical reaction networks. Journal of Mathematical Chemistry, 50:274–288, 2012.
21. M. D. Johnston, D. Siegel, and G. Szederkényi. Computing weakly reversible linearly conjugate chemical reaction networks with minimal deficiency. Mathematical Biosciences, 241:88–98, 2013.
22. W. Liebermeister and E. Klipp. Biochemical networks with uncertain parameters. IEE Proceedings Systems Biology, 152:97–107, 2005.
23. G. Lipták, G. Szederkényi, and K. M. Hangos. Computing zero deficiency realizations of kinetic systems. Systems & Control Letters, 81:24–30, 2015.
24. G. Lipták, G. Szederkényi, and K. M. Hangos. Kinetic feedback design for polynomial systems. Journal of Process Control, 41:56–66, 2016.
25. L. Ljung. System Identification - Theory for the User. Prentice Hall, 1999.
26. F. Llaneras and J. Picó. An interval approach for dealing with flux distributions and elementary modes activity patterns. Journal of Theoretical Biology, 246:290–308, 2007.
27. H. F. Lodish. Molecular cell biology. W.H. Freeman, New York, 2000.
28. T. Nagy and T. Turányi. Uncertainty of Arrhenius parameters. International Journal of Chemical Kinetics, 43:359–378, 2011.
29. C. Schillings, M. Sunnaker, J. Stelling, and C. Schwab. Efficient characterization of parametric uncertainty of complex biochemical networks. PLOS Computational Biology, 11(8):e1004457, 2015.
30. S. Schnell, M. J. Chappell, N. D. Evans, and M. R. Roussel. The mechanism distinguishability problem in biochemical kinetics: The single-enzyme, single-substrate reaction as a case study. Comptes Rendus Biologies, 329:51–61, 2006.
31. G. Shinar and M. Feinberg. Structural sources of robustness in biochemical reaction networks. Science, 327:1389–1391, 2010.
32. E. Sontag. Structure and stability of certain chemical networks and applications to the kinetic proofreading model of T-cell receptor signal transduction. IEEE Transactions on Automatic Control, 46:1028–1047, 2001.
33. G. Szederkényi. Computing sparse and dense realizations of reaction kinetic systems. Journal of Mathematical Chemistry, 47:551–568, 2010.
34. G. Szederkényi, J. R. Banga, and A. A. Alonso. Inference of complex biological networks: distinguishability issues and optimization-based solutions. BMC Systems Biology, 5:177, 2011.
35. G. Szederkényi and K. M. Hangos. Finding complex balanced and detailed balanced realizations of chemical reaction networks. Journal of Mathematical Chemistry, 49:1163–1179, 2011.
36. G. Szederkényi, K. M. Hangos, and T. Péni. Maximal and minimal realizations of reaction kinetic systems: computation and properties. MATCH Commun. Math. Comput. Chem., 65:309–332, 2011.
37. G. Szederkényi, B. Ács, and G. Szlobodnyik. Structural analysis of kinetic systems with uncertain parameters. In 2nd IFAC Workshop on Thermodynamic Foundations for a Mathematical Systems Theory - TFMST, IFAC Papersonline, volume 49, pages 24–27, 2016.
38. Y. Takeuchi. Global Dynamical Properties of Lotka-Volterra Systems. World Scientific, Singapore, 1996.
39. Z. A. Tuza and G. Szederkényi. Computing core reactions of uncertain polynomial kinetic systems. In 23rd Mediterranean Conference on Control and Automation (MED), June 16-19, 2015. Torremolinos, Spain, pages 1187–1194, 2015.
40. A. Weinmann. Uncertain Models and Robust Control. Springer-Verlag, 1991.
41. J-L. Wu. Robust stabilization for single-input polytopic nonlinear systems. IEEE Transactions on Automatic Control, 2006.
42. T. M. Yi, H. Kitano, and M. I. Simon. A quantitative characterization of the yeast heterotrimeric g protein cycle. Proc. Natl. Acad. Sci. USA, 100(19):10764–10769, 2003.
43. B. Ács, G. Szederkényi, and D. Csercsik. A new efficient algorithm for determining all structurally different realizations of kinetic systems. MATCH Commun. Math. Comput. Chem., 77:299–320, 2017.
44. B. Ács, G. Szederkényi, Z. A. Tuza, and Z. Tuza. Computing linearly conjugate weakly reversible kinetic structures using optimization and graph theory. MATCH Commun. Math. Comput. Chem., 74:481–504, 2015.
45. B. Ács, G. Szederkényi, Zs. Tuza, and Z. A. Tuza. Computing all possible graph structures describing linearly conjugate realizations of kinetic systems. Computer Physics Communications, 204:11–20, 2016.
46. P. Érdi and J. Tóth. Mathematical Models of Chemical Reactions. Theory and Applications of Deterministic and Stochastic Models. Manchester University Press, Princeton University Press, Manchester, Princeton, 1989.