Inference of complex biological networks: distinguishability issues and optimizationbased solutions
 Gábor Szederkényi^{1, 2}Email author,
 Julio R Banga^{1}Email author and
 Antonio A Alonso^{1}
DOI: 10.1186/175205095177
© Szederkényi et al; licensee BioMed Central Ltd. 2011
Received: 31 May 2011
Accepted: 28 October 2011
Published: 28 October 2011
Abstract
Background
The inference of biological networks from highthroughput data has received huge attention during the last decade and can be considered an important problem class in systems biology. However, it has been recognized that reliable network inference remains an unsolved problem. Most authors have identified lack of data and deficiencies in the inference algorithms as the main reasons for this situation.
Results
We claim that another major difficulty for solving these inference problems is the frequent lack of uniqueness of many of these networks, especially when prior assumptions have not been taken properly into account. Our contributions aid the distinguishability analysis of chemical reaction network (CRN) models with mass action dynamics. The novel methods are based on linear programming (LP), therefore they allow the efficient analysis of CRNs containing several hundred complexes and reactions. Using these new tools and also previously published ones to obtain the network structure of biological systems from the literature, we find that, often, a unique topology cannot be determined, even if the structure of the corresponding mathematical model is assumed to be known and all dynamical variables are measurable. In other words, certain mechanisms may remain undetected (or they are falsely detected) while the inferred model is fully consistent with the measured data. It is also shown that sparsity enforcing approaches for determining 'true' reaction structures are generally not enough without additional prior information.
Conclusions
The inference of biological networks can be an extremely challenging problem even in the utopian case of perfect experimental information. Unfortunately, the practical situation is often more complex than that, since the measurements are typically incomplete, noisy and sometimes dynamically not rich enough, introducing further obstacles to the structure/parameter estimation process. In this paper, we show how the structural uniqueness and identifiability of the models can be guaranteed by carefully adding extra constraints, and that these important properties can be checked through appropriate computation methods.
Background
During the last decade, the wide availability of highthroughput biological data has made it possible to produce new knowledge via a systems biology approach [1–3]. The inference of biochemical networks (i.e. the mathematical mapping of the molecular interactions in the cell) is therefore a question of key importance in the field. During the last decade, many methods have been developed to solve the networkinference (sometimes called reverseengineering [4]) problems arising in e.g. gene expression [5–13], signal transduction [14–17] and metabolic networks [18–25].
In this context, it is particularly worth mentioning the DREAM initiative (Dialogue for Reverse Engineering Assessments and Methods) [26], which targeted the problems of cellular network inference and quantitative model building in systems biology. DREAM tries to address two fundamental questions: (i) how can we assess how well we are describing the networks of interacting molecules that underlie biological systems? and (ii) how can we know how well we are predicting the outcome of previously unseen experiments from our models? Interestingly, one of the main conclusions of the DREAM3 event was that the vast majority of the teams' predictions were statistically equivalent to random guesses. Moreover, even for particular problem instances like gene regulation network inference, there was no onesizefitsall algorithm [27].
The use of a performance profiling framework with the DREAM3 benchmark problems revealed that current inference methods are affected by different types of systematic prediction errors [6]. These authors conclude that reliable network inference from gene expression data remains an unsolved problem. Further, they highlight two major difficulties in the case of genenetwork reverse engineering: limited data (which may leave the inference problem underdetermined), and the difficulty of distinguishing direct from indirect regulation. Prill et al [27] further explored the issue of intrinsic impediments to network inference, designating identifiability of certain network edges and systematic false positives as the main barriers. In this paper, we consider the widely used reaction kinetic formalism, where dynamic models of biological networks are described by a set of ordinary differential equations (see, e.g. [28–30] and the related literature). In particular, we consider the central question of the identifiability of such a network as understood in the systems and control area [31, 32].
Identifiability analysis studies whether there is a theoretical chance of uniquely determining the parameters of a mathematical model assuming perfect noisefree measurements and errorfree modeling [33–35]. One of the early approaches for identifiability testing of nonlinear models is based on the Taylorseries expansion of the system output using the fact that the Taylor coefficients are unique [36]. A similar but more general method uses the generating series or Volterraseries coefficients of the system which is the nonlinear generalization of the Laplacetransform method used for linear systems [37]. In [38] a similarity transformation approach is proposed that gives necessary and sufficient conditions on local and global identifiability through the checking of nonlinear controllability and observability conditions. The appearance of differential algebra methods in systems and control theory [39, 40] opened the possibility for new types of identifiability tests that have gained significant popularity [41–43]. Further theoretical developments in the field include the identifiability conditions of rational function state space models [43], the possible effect of initial conditions on identifiability [44], and the application of Liealgebras [45]. While identifiability is the property of a certain parameterized model, a related notion called distinguishability addresses the problem whether two or more parameterized models (with the same or with different structure) can produce the same output for any allowed input [46–48]. The literature about identifiability and distinguishability of biological and chemical system models is relatively wide: Compartmental systems (that form a special subclass of general massaction networks) are studied in [38, 49, 50]. The authors treat general nonlinear CRNs in [51, 52] and [53] where it is shown that for thermodynamically meaningful models, nonlinearity reduces the chance of indistinguishability compared to the linear case [54]. Geometric conditions for the indistinguishability of CRNs are given in [55] with a related comment in [56]. Computer algebra tools can be successfully used for the symbolic computations needed for identifiability and distinguishability testing of complex models [57–60].
The importance of identifiability has been recognized previously in systems biology, too [14, 61–64]. However, and despite a number of works illustrating ways to test the structural and practical identifiability of models [65–67], a significant portion of modeling studies in systems biology continue to ignore this key property.
It has been known for long that chemical reaction networks with different structure and/or parametrization may produce the same dynamical models describing the timeevolution of species concentrations [28, 55]. A related problem, namely the nonunique structure of Petri nets associated to reaction network dynamics, is studied in [68]. Additionally, the value of prior information in biological network inference was clearly shown in [69, 70] by applying Bayesian network models. However, a constructive optimizationbased approach for the study of dynamically equivalent (or similar) reaction networks is a recent development [71–74], which we further extend in this paper.
As a novelty, we present in this paper the definition and a computational method to find the socalled core reactions that are present in any dynamically equivalent reaction network if the set of complexes is given a priori. Moreover, a computationally improved method is introduced for the computation of dense realizations of CRNs together with a modified algorithm to check the uniqueness of a constrained reaction network structure. Structural nonuniqueness and the use of the proposed computational methods will be illustrated with the help of biological models known from the literature.
The structure of the paper is the following. The 'Methods' section introduces the notions of chemical reaction networks, structural identifiability and distinguishability of dynamical models. Moreover, it contains the procedures to obtain core reactions of a network and its sparse and dense representations, which rely on standard methods of linear programming (LP) and mixed integer linear programming (MILP) [75–78]. The analysis of four biological system models can be found in the 'Results and discussion' section, followed by the conclusions.
Methods
where x ∈ ℝ ^{ n } is the state vector, y ∈ ℝ ^{ m } is the output, u ∈ ℝ ^{ k } is the input, and θ ∈ ℝ ^{ d } denotes the parameter vector. We assume that the functions f and h are polynomial in the variables x, u and θ. Clearly, mass action type CRNs described in the following subsection (where θ is typically the set of reaction rate coefficients), and simple deterministic models of gene regulation such as the one in Example 4 belong to this model class.
Basic notions and known results related to massaction models
In this subsection, the basic definitions for the description of CRNs will be given together with the already published results on finding dynamically equivalent network realizations with certain prescribed properties.
Structural and dynamical description of massaction networks
 1.
$\mathcal{S}=\left\{{X}_{1},\dots ,{X}_{n}\right\}$ is the set of species or chemical substances.
 2.$\mathcal{C}=\left\{{C}_{1},\dots ,{C}_{m}\right\}$ is the set of complexes. Formally, the complexes are represented as linear combinations of the species, i.e.${C}_{i}=\sum _{j=1}^{n}{\alpha}_{ij}{X}_{j},\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}i=1,\dots ,m,$(2)
where α_{ ij } are nonnegative integers and are called the stoichiometric coefficients.
 3.
$\mathcal{R}=\{\left({C}_{i},{C}_{j}\right){C}_{i},{C}_{j}\in \mathcal{C}$, and C_{ i } is transformed to C_{ j } in the CRN} is the set of reactions. The relation $\left({C}_{i},{C}_{j}\right)\in \mathcal{R}$ will be denoted as C_{ i } → C_{ j } . Moreover, a nonnegative weight, the reaction rate coefficient denoted by k_{ ij } is assigned to each reaction C_{ i } → C_{ j } . Naturally, if the reaction C_{ i } → C_{ j } is not present in the CRN then k_{ ij } = 0.
The above characterization naturally gives rise to the following graph structure (often called 'FeinbergHornJackson graph' or simply reaction graph) of a CRN [29]. The weighted directed graph G = (V, E) of a CRN consists of a finite nonempty set V of vertices and a finite set E of ordered pairs of distinct vertices called directed edges. The vertices correspond to the complexes, i.e. V = {C_{1}, C_{2}, ... C_{ m } }, while the directed edges represent the reactions, i.e. (C_{ i } , C_{ j } ) ∈ E if complex C_{ i } is transformed to C_{ j } in the CRN. The positive reaction rate coefficients k_{ ij } are assigned as weights to the corresponding directed edges C_{ i } → C_{ j } in the graph. (Edges corresponding to zero rate coefficients are not drawn in the reaction graph.) A set of complexes {C_{1}, ..., C_{ k } } is called a linkage class of a CRN, if the complexes of the set are linked to each other in the reaction graph but not to any other complex. It is remarked that loops (i.e. directed edges that start and end at the same vertex) are not allowed in reaction graphs.
Dynamical equivalence of massaction networks
In this case, $\left({Y}^{\left(i\right)}{A}_{k}^{\left(i\right)}\right)$ for i = 1, 2 are called realizations of a kinetic vector field f (see, e.g. [80] for more details). It is also appropriate to call $\left({Y}^{\left(1\right)},{A}_{k}^{\left(1\right)}\right)$ a realization of $\left({Y}^{\left(2\right)},{A}_{k}^{\left(2\right)}\right)$ and vice versa.
where ${A}_{k}^{\left(1\right)}$ and ${A}_{k}^{\left(2\right)}$ are valid Kirchhoff matrices and M is the invariant matrix containing the coefficients of the monomials.
Among the dynamically equivalent realizations, it is important to recall the following characteristic ones described in [71, 72]. A sparse realization contains the minimal number of reactions that is needed for the exact description of the corresponding dynamics (3). A dense realization contains the maximal number of reactions among dynamically equivalent realizations with a fixed complex set given by Y. While sparse realizations are generally structurally nonunique (as it will be illustrated for the constrained case, too, in Example 1), the structure of dense realizations with a given complex set is unique, and it contains every possible dynamically equivalent structure as a proper subgraph (i.e. a dense realization is a kind of superstructure) [71].
Known computation approaches for finding preferred CRN realizations
Here we briefly summarize the already published results corresponding to the computation of preferred dynamically equivalent CRN realizations (more details of these methods can be found in the publications [71–73, 81]). The computation of dense and sparse realizations can be traced back to mixed integer linear programming (MILP) where the decision variables are the nondiagonal elements of A_{ k } , the linear constraints encode the kinetic properties of the model, and the objective function contains integer variables for minimizing/maximizing the number of nonzero reaction rate coefficients [72]. It is remarked that the computation of sparse realizations is an NPhard problem, where generally mixed integer linear programming cannot be avoided [82]. There exist certain conditions under which the problem can be solved in polynomial time [83] but these are often not fulfilled in the case of CRNs. Moreover, there are effective heuristics to address the problem [84], but convergence to one of the truly sparsest structures is not guaranteed. Luckily, the MILPbased computation of sparse CRN realizations can be parallelized effectively thus allowing a larger number of complexes to be treated. The computation of realizations having the minimal/maximal number of complexes or the reversibility property can also be solved in the MILP framework [71]. Moreover, it was shown in [73] that finding detailed balanced and complex balanced realizations of CRNs is a simple linear programming (LP) problem. Finally, weakly reversible dynamically equivalent CRN realizations can also be determined (if they exist) using MILP [85].
Constrained realizations of CRNs and testing their structural uniqueness

P1 Given a CRN (Y, A_{ k } ) and a constraint set $\mathcal{K}$, the unweighted reaction graph of any constrained realization is the subgraph of the unweighted reaction graph of the constrained dense realization.

P2 If the sets of complexes and constraints are fixed, then for any CRN, the structure of the constrained dense realization is unique.

P3 The reaction graph structure of a CRN with given sets of complexes and constraints is unique if and only if the unweighted directed graphs of its constrained dense and sparse realizations are identical.
The proofs of P1, P2 and P3 follow similar (although not completely identical) lines that were published in [71], and they are given for convenience in the Appendix at the end of the paper.
New concepts and computation results related to dynamically equivalent networks
This subsection contains new methodological contributions that extend the previously published results.
Making the computation of dense realizations more efficient
Computing dense realizations is treated originally also in a MILPframework in [72]. However, using the structural uniqueness of such realizations given by P1, it is easy to give a polynomialtime algorithm based on a finite series of linear programming (LP) optimization steps. The idea of the improved algorithm is simple: the reaction C_{ i } → C_{ j } belongs to the (constrained) dense realization if and only if there exists any dynamically equivalent (constrained) realization where [A_{ k } ] _{ j,i } > 0. The result directly follows from the fact that the unweighted reaction graphs of (constrained) dense realizations give a unique superstructure. This allows us to formulate a polynomialtime method based on pure LP to determine (constrained) dense realizations as follows.
It is important to remark that the definition of ε_{ ij } in the form of (11) guarantees the solvability of (12). Naturally, the above described method can also be used for determining constrained dense realizations by adding constraints of the form (9) to the LP problems (10) and (12).
Using the notion and described properties of constrained realizations, we are now able to test the structural uniqueness of given CRNs. To accomplish this, only the (constrained) dense and sparse realizations have to be computed and compared (see P3). This method will be illustrated in Example 2.
Definition and computation of core and noncore reactions
where the matrix A_{ k } contains the decision variables, and the known matrices are Y and M. It is wellknown that this task is equivalent to an LP problem where the objective function is an arbitrary linear function of the elements of A_{ k } [76]. Then, reaction C_{ p } → C_{ q } is a corereaction if and only if the set defined by (13)(17) is empty (i.e. the corresponding LP problem is infeasible), because in this case there is no dynamically equivalent CRN realization where C_{ p } → C_{ q } is not present. We remark here that the presented procedures for determining constrained dense realizations and computing core reactions are parallel in their original forms since the individual LP steps are independent of each other. Therefore the proposed methods can be very effectively implemented in a grid or multicore hardware environment [86].
Basic concepts on structural identifiability and distinguishability
and x(t, θ) denotes the solution of (1) with parameter vector θ. According to (18), a structurally nonidentifiable model can produce exactly the same observed output with different parametrization. This is clearly a fundamental obstacle of determining the true model parameters from measurements even if the selected model structure is considered to be correct.
Let us denote two parameterized models with possibly different structure by ${\mathcal{M}}_{1}\left({\theta}_{1}\right)$ and ${\mathcal{M}}_{2}\left({\theta}_{2}\right)$, respectively, where θ_{ i } denote the parameter vector. Then ${\mathcal{M}}_{1}$ is called distinguishable from ${\mathcal{M}}_{2}$ if for any θ_{1} (possibly except for a finite number of values) there is no θ_{2} such that the inputoutput behaviour of ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$ is the same [47]. Clearly, if ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$ are indistinguishable and both model structures are feasible in a certain application, then there is no way to decide from inputoutput measurements to which one corresponds to the true model that generated the data.
In the case of CRNs, we will assume that all species concentrations are measured (i.e. y = x), the input is zero (i.e. we study autonomous systems), and that the set of possible chemical complexes is given. Generally, the model parameter vector θ is the set of reaction rate coefficients which are the offdiagonal elements of A_{ k } . Clearly, if a CRN has several different dynamically equivalent realizations, then these realizations are not distinguishable without additional constraints, and the model cannot be identifiable if all the rate coefficients are to be determined [55]. This situation can be improved by using prior knowledge in the form of adding further constraints on the model parameters such as the simple ones given by eq. (9). This way, the number of parameters to be estimated can be reduced and/or their feasibility region can be shrinked. It is important to note that although the structural uniqueness of a CRN definitely reduces the degree of nonidentifiability (since zero and nonzero parameters are separated), it does not necessarily imply structural identifiability [53], and this latter property has to be checked by further numerical methods [31, 32].
Results and discussion
In this section, the application of the previously mentioned methods for finding different dynamically equivalent structures will be illustrated using biological models taken from the literature. The detailed numerical data corresponding to Examples 13 are contained in a standard spreadsheet form with brief explanations in Additional file 1: CRN_data.xls.
Example 1: a positive feedback motif
Our starting point is that we have a dynamic model of the process in the standard polynomial form of (20)(24), the parameters of which are known from the results of identification and/or from literature. As we will see below, without welldefined constraints on the possible set of complexes and reactions, exactly the same dynamics can be realized in principle by a wide range of mechanisms.
It is visible that the structure of the dense realization is the same as in Figure 1b but the core reactions are different from the ones shown there. Here the degradation of mRNA is described by a core reaction but interestingly, the reaction corresponding to translation is not a core one. Naturally, this implies that the possible sparse realization structures with the second parametrization are different from the ones shown in Figure 2. Note that here the only goal was to illustrate the possible change of core and noncore reactions, and therefore the biological relevance of the parameter values in eq. (28) is not assumed in this case.
The dense CRN realization of the dynamics (20)(24) with the updated Y' matrix given in eq. (29) using the original parameters described in (27) is shown in Figure 3b, where the core and noncore reactions are again indicated. It is apparent that now there are only 5 core reactions, and none of the remaining 12 reactions are essential to represent the dynamics (20)(24). This means that the introduction of a new complex increased the flexibility of the network (i.e. mathematically, the majority of the reactions can be substituted by other ones and the network still maintains its original dynamics). Of course, not any combination of the noncore reactions can be omitted from the network, because the sparse realizations show that at least 9 reactions are needed to keep dynamical equivalence. It can be computed easily that the theoretical maximum number of sparse realizations with different structures is $\left(\begin{array}{c}\hfill 12\hfill \\ \hfill 179\hfill \end{array}\right)=495$. However, as the numerical experiments show, majority of these structures do not give a practically feasible dynamically equivalent realization.
The above results clearly show that certain mechanisms may remain undetectable (or they are falsely detected) even if we have complete species concentration measurements and full information about possible complex formation, that are not very realistic assumptions. Moreover, the sparsest dynamically equivalent structure of massaction models is not unique, therefore sparsity enforcing approaches for determining 'true' reaction structures are not enough in themselves without the necessary amount of prior information given in the form of additional constraints. The practical situation is most often even worse than that, since the measurements are typically incomplete, noisy and sometimes dynamically not rich enough, that may introduce further obstacles to the structure/parameter estimation process [66, 89].
Example 2: a biochemical switch in yeast cells
The following example is taken from [90] and it describes a 'switching device' in yeast cycle regulation. The detailed system description can be found in [90] and in the accompanying supporting information document. The order of state variables, corresponding to concentrations, is the same as in the original article, and is shown below:
The structure of the dense realization indicating the 12 core and 16 noncore reactions can be seen in Figure 4b.
It can be shown using the computational methods described in the 'Methods' section that the only possible sparse realization structure is identical to that of the original network. Therefore in this special case, there is only one possible reaction structure containing the minimal number of reactions. A straightforward approach to ensure the structural uniqueness of the whole network is to exclude all reactions that are not meaningful from the examined application's point of view or that are contradictory to modeling assumptions. For the current example, the removal of an unexpectedly low number of reactions is enough to obtain a unique structure. It can be shown by computing the corresponding constrained dense and sparse realizations, that excluding the reactions X_{5} → X_{3} + X_{5}, X_{4} → X_{3} + X_{4}, X_{2} + X_{3} → X_{3} + X_{5}, and X_{3} + X_{1} → X_{3} + X_{4} is enough to make the reaction structure unique that is identical to the original structure shown in Figure 4. In other words, the exclusion of 4 wellselected reactions leads to the removal of an additional 6 reactions leaving only 18.
Example 3: a repressilator structure with 5 nodes and autoactivation
for the index pairs (i, j) ∈ {(1, 5), (2, 1), (3, 2), (4, 3), (5, 4)}. In eqs. (33)(39), G_{ i } and P_{ i } represent the ith gene and protein, respectively. For the genes, superscripts A and R refer to activated and repressed states, respectively. Let us denote with r_{ i,k } the reaction with rate coefficient k_{ i,k } in eqs. (33)(39).
In particular, it is easy to show (see also Additional File 1: CRN_data.xls) that reactions ${G}_{i}^{AR}\to {G}_{i}^{R}+{P}_{i}$ and ${G}_{i}^{AR}\to {G}_{i}^{R}$ are always indistinguishable. Similarly, the reaction ${G}_{i}^{R}+{P}_{i}\to {G}_{i}^{AR}$ can be substituted with the combination of reactions ${G}_{i}^{R}+{P}_{i}\to {G}_{i}^{R}$ and ${G}_{i}^{R}+{P}_{i}\to {G}_{i}^{AR}+{P}_{i}$. It can be seen from these results that in order to have a model with unique structure, it is very important to a priori exclude all reactions that are not meaningful for the particular application.
Example 4: sparse linear gene regulation network models
where A ∈ ℝ ^{ n×n } contains the connectivity information of the network. A_{ i,j } > 0 indicates activation from node j to node i, while A_{ i,j } < 0 means repression, diagonal elements of A represent autoactivation or autorepression depending on their sign. x ∈ ℝ ^{ n } is the fully or partially measurable state of the system describing the time evolution of concentrations, and the input part Bu represents experimental perturbation (e.g. activation) of the genes. It is also a common assumption that the network is 'sparse' which means that there are only a limited number of activation or repression links between the nodes (i.e. the matrix A is 'sparse', too). But assuming sparsity can be a serious obstacle to identifiability as it will be shown.
where '+', '' and '*' represent positive, negative and nonzero (but otherwise undefined) parameter values, respectively. If there are no prior assumptions about the structure of the interconnection matrix or about the relations between certain parameters, we can easily test the structural nonidentifiability of the model by checking whether all nodes are reachable from the perturbed node on a directed path in the interconnection graph or not [33]. The reachability of nodes can be tested by several methods, e.g. a depth or breadthfirstsearch (DFS or BFS) of the corresponding directed graphs that are fast polynomialtime algorithms [94]. To give a very simple example, it is clear from Figure 6a, that if nodes 1 or 2 were excited by an input signal, then the connections between the other nodes (310) would be undetectable by any method, however sophisticated it is. To examine whether situations like this one are common, we generated 10000 random state space models using the same method, and assuming zero initial conditions as in [93]. The connectivity of the corresponding directed graphs was tested using DFS. For 10 nodes and 2 nonzero elements in each row of A (i.e. N = 10, K = 2), we obtained that 73.38% of the generated models are structurally nonidentifiable. The histogram showing the number of reachable states is shown in Figure 6b. The situation is dramatically improving if K is increased to 3 as shown in Figure 6c. In this case, around 17% of the models are structurally nonidentifiable. When we have 20 nodes and 5 nonzero elements in each row of A (the second case investigated in [93]), then only 1.6% of the generated models are structurally nonidentifiable. The results show that 'sparsity' has a clearly negative effect on structural identifiability because of limited information transmission between nodes. And finally, we did not speak at all about practical identifiability which is known to be a challenging issue even if the required structural properties are fulfilled [66].
Relation between high level networks and CRN structure
As shown in Example 3, the various possible dynamically equivalent CRN structures do not correspond to a different GRN structure, if all species concentration measurements are available and the mapping described in [92] is used for transforming the models into each other. Hence, exact matching of the dynamics of different GRN structures may generally be a too severe restriction. To extend this line of research, the relaxation of dynamical equivalence to 'close dynamical similarity' seems to be more meaningful but the corresponding definitions and computational methods are much more complex than in the case of dynamical equivalence. One promising recent approach to assess dynamical similarity of CRNs (that also adds more degrees of freedom to the computations) is the concept of 'linear conjugacy' [74]. However, it might happen that dynamically completely equivalent GRN structures will be shown in the future.
Conclusions
It has been shown in this paper using illustrative examples that biological network structures modeled by CRNs often cannot be uniquely determined even if the structure of the corresponding mathematical model is assumed to be known and all dynamical variables are measurable. The structural uniqueness and identifiability of the models often require additional constraints.
The main new contributions of the paper are the following. Firstly, core reactions present in any dynamically equivalent CRN realizations with a given complex set have been defined and a simple procedure with polynomial timecomplexity has been given for determining them. Clearly, the core reactions are mandatory elements of every dynamically equivalent CRN realization assuming a fixed complex set. Secondly, a polynomialtime method based on linear programming for computing dense realizations has been outlined that is more scalable and therefore presents a clear improvement over the previously used MILPbased method. As an additional minor extension of previous results, constrained realizations of CRNs have been defined, and a computational method has been proposed to check the uniqueness of constrained realizations.
The presented concepts and algorithms were illustrated on previously published models describing biological processes. It was shown that the set of core reactions may change with the modification of the complex set. The examples also show that the frequently applied sparsity assumption alone is not enough for structural uniqueness of CRNs. Moreover, in the case of simple linear genetic network models, too sparse structures can degrade identifiability properties. The results further support the fact that as much prior information as possible should be incorporated in structural and parametric inference problems.
A Appendix
where $\u0232\in {\mathbb{R}}^{\left(n+z\right)\times m}$ and $\stackrel{\u0304}{b}\in {\mathbb{R}}^{n+z}$. Let us assume now that p is a dense solution for (48), i.e. it contains the maximal possible number of nonzero elements. If p has no zero elements, then the result to be proved is trivially satisfied. Therefore, without the loss of generality we can assume that the first l < m elements of p are nonzero, while the rest are zero, i.e. p_{ j } ≠ 0 for j = 1, ..., l, and p_{ j } = 0 for j = l + 1, ..., m. This can always be achieved by the appropriate reordering of the elements of p. Assume now that p' ∈ ℝ ^{ m } is also a solution for (48), but ${p}_{c}^{\prime}\ne 0$ for some c ∈ ℤ, l + 1 ≤ c ≤ m. Then p' = p + v, where $\u0232\cdot v=0$, and v_{ c } ≠ 0. In this case, p″ = p + λ · v is also a solution for (48) for any λ ∈ ℝ and λ can always be chosen so that ${p}_{j}^{\u2033}\ne 0$ for j = 1, ..., l, and there is at least one index l + 1 ≤ c ≤ m for which ${p}_{c}^{\u2033}\ne 0$. However, this contradicts to the assumption that p is a dense solution for (48).
Proof of P2. This is a straightforward consequence of P1, since the unweighted directed graphs of all constrained dense realizations must be identical.
Proof of P3. If the graph structure of the constrained realization is unique, then it trivially implies that the structures of the constrained dense and sparse realizations are identical, since there exists only one possible constrained reaction structure. If the structures of the constrained dense and sparse realizations are identical, then the number of nonzero reaction rates is the same in any constrained realizations including the constrained dense ones. Then it follows from P1 that the constrained reaction structure is unique.
Declarations
Acknowledgements
This work was financially supported by project CAFE (Computer Aided Food Process Engineering) FP7KBBE20071 (Grant no: 212754). The authors acknowledge the support by the Spanish government, MICINN project 'MultiSysBio' (ref. DPI200806880C0302), and by CSIC intramural project 'BioREDES' (ref. PIE201170E018). GS acknowledges the partial support of the Hungarian Scientific Research Fund through grant no. K83440. The authors thank Dr. Irene Otero Muras (Dept. of Computational Systems Biology, ETH Zurich) and Prof. Zsolt Tuza (MTA SZTAKI) for the fruitful discussions. The authors wish to thank the anonymous reviewers for their helpful comments.
Authors’ Affiliations
References
 Wolkenhauer O: Systems biology: The reincarnation of systems theory applied in biology?. Briefings in Bioinformatics. 2001, 2: 258270. 10.1093/bib/2.3.258.View ArticlePubMed
 Stelling J: Mathematical models in microbial systems biology. Current Opinion in Microbiology. 2004, 7: 513518. 10.1016/j.mib.2004.08.004.View ArticlePubMed
 Kitano H: Computational systems biology. Nature. 2002, 420: 206210. 10.1038/nature01254.View ArticlePubMed
 Csete M, Doyle J: Reverse engineering of biological complexity. Science. 2002, 295: 16641669. 10.1126/science.1069981.View ArticlePubMed
 De Jong H: Modeling and simulation of genetic regulatory systems: a literature review. Journal of Computational Biology. 2002, 9: 67103. 10.1089/10665270252833208.View ArticlePubMed
 Marbach D, Prill R, Schaffter T, Mattiussi C, Floreano D, Stolovitzky G: Revealing strengths and weaknesses of methods for gene network inference. Proceedings of the National Academy of Sciences of the United States of America. 2010, 107: 62866291. 10.1073/pnas.0913357107.PubMed CentralView ArticlePubMed
 Ay A, Arnosti D: Mathematical modeling of gene expression: a guide for the perplexed biologist. Critical Reviews in Biochemistry and Molecular Biology. 2011, 46: 137151. 10.3109/10409238.2011.556597.PubMed CentralView ArticlePubMed
 Hecker M, Lambeck S, Toepfer S, Van Someren E, Guthke R: Gene regulatory network inference: data integration in dynamic models  a review. Biosystems. 2009, 96: 86103. 10.1016/j.biosystems.2008.12.004.View ArticlePubMed
 Bansal M, Belcastro V, AmbesiImpiombato A, Di Bernardo D: How to infer gene networks from expression profiles. Molecular Systems Biology. 2007, 3: article number 78
 Tegner J, Yeung M, Hasty J, Collins J: Reverse engineering gene networks: integrating genetic perturbations with dynamical modeling. Proceedings of the National Academy of Sciences of the United States of America. 2003, 100: 59445949. 10.1073/pnas.0933416100.PubMed CentralView ArticlePubMed
 de la Fuente A, Brazhnik P, Mendes P: Linking the genes: inferring quantitative gene networks from microarray data. Trends in Genetics. 2002, 18: 395398. 10.1016/S01689525(02)026926.View ArticlePubMed
 Ronen M, Rosenberg R, Shraiman B, Alon U: Assigning numbers to the arrows: parameterizing a gene regulation network by using accurate expression kinetics. Proceedings of the National Academy of Sciences of the United States of America. 2002, 99: 1055510560. 10.1073/pnas.152046799.PubMed CentralView ArticlePubMed
 Thomas R, Paredes CJ, Mehrotra S, Hatzimanikatis V, Papoutsakis ET: A modelbased optimization framework for the inference of regulatory interactions using timecourse DNA microarray expression data. BMC Bioinformatics. 2007, 8: 228239. 10.1186/147121058228.PubMed CentralView ArticlePubMed
 Schaber J, Klipp E: Modelbased inference of biochemical parameters and dynamic properties of microbial signal transduction networks. Current Opinion in Biotechnology. 2010, 22: 109116.View ArticlePubMed
 Li S, Assmann S, Albert R: Predicting essential components of signal transduction networks: a dynamic model of guard cell abscisic acid signaling. PLoS Biology. 2006, 4: 17321748.View Article
 SaezRodriguez J, Kremling A, Conzelmann H, Bettenbrock K, Gilles E: Modular analysis of signal transduction networks. IEEE Control Systems Magazine. 2004, 24: 3552.View Article
 Klamt S, SaezRodriguez J, Lindquist J, Simeoni L, Gilles E: A methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinformatics. 2006, 7: 5610.1186/14712105756. (pp. 126)PubMed CentralView ArticlePubMed
 Arkin A, Shen P, Ross J: A test case of correlation metric construction of a reaction pathway from measurements. Science. 1997, 277: 12751279. 10.1126/science.277.5330.1275.View Article
 Stelling J, Klamt S, Bettenbrock K, Schuster S, Gilles E: Metabolic network structure determines key aspects of functionality and regulation. Nature. 2002, 420: 190193. 10.1038/nature01166.View ArticlePubMed
 Jeong H, Tombor B, Albert R, Oltvai Z, Barabási A: The largescale organization of metabolic networks. Nature. 2000, 407: 651654. 10.1038/35036627.View ArticlePubMed
 Forster J, Famili I, Fu P, Palsson B, Nielsen J: Genomescale reconstruction of the Saccharomyces cerevisiae metabolic network. Genome Research. 2003, 13: 24410.1101/gr.234503.PubMed CentralView ArticlePubMed
 Kell D: Metabolomics and systems biology: making sense of the soup. Current Opinion in Microbiology. 2004, 7: 296307. 10.1016/j.mib.2004.04.012.View ArticlePubMed
 Sauer U: Metabolic networks in motion: 13Cbased flux analysis. Molecular Systems Biology. 2006, 2: article number 62
 Crampin E, Schnell S, McSharry P: Mathematical and computational techniques to deduce complex biochemical reaction mechanisms. Progress in Biophysics and Molecular Biology. 2004, 86: 77112. 10.1016/j.pbiomolbio.2004.04.002.View ArticlePubMed
 Yuan Y, Stan GB, Warnick S, Goncalves J: Robust dynamical network structure reconstruction. Automatica. 2011, 47: 12301235. 10.1016/j.automatica.2011.03.008.View Article
 Stolovitzky G, Monroe D, Califano A: Dialogue on ReverseEngineering Assessment and Methods. Annals of the New York Academy of Sciences. 2007, 1115: 122. 10.1196/annals.1407.021.View ArticlePubMed
 Prill R, Marbach D, SaezRodriguez J, Sorger P, Alexopoulos L, Xue X, Clarke N, AltanBonnet G, Stolovitzky G: Towards a rigorous assessment of systems biology models: the DREAM3 challenges. PLoS ONE. 2010, 5: e920210.1371/journal.pone.0009202. (pp. 118)PubMed CentralView ArticlePubMed
 Horn F, Jackson R: General mass action kinetics. Archive for Rational Mechanics and Analysis. 1972, 47: 81116.View Article
 Feinberg M: Chemical reaction network structure and the stability of complex isothermal reactors  I. The deficiency zero and deficiency one theorems. Chemical Engineering Science. 1987, 42: 22292268. 10.1016/00092509(87)800994.View Article
 Érdi P, Tóth J: Mathematical Models of Chemical Reactions. Theory and Applications of Deterministic and Stochastic Models. 1989, Manchester, Princeton: Manchester University Press, Princeton University Press
 Ljung L: System Identification: Theory for the User. 1999, Prentice Hall, 2
 Walter E, Pronzato L: Identification of Parametric Models. 1997, Springer
 Bellmann R, Aström KJ: On structural identifiability. Mathematical Biosciences. 1970, 7: 329339. 10.1016/00255564(70)90132X.View Article
 Walter E: Identification of State Space Models. 1982, SpringerView Article
 Walter E: Identifiability of Parametric models. 1987, Pergamon Press, Oxford
 Pohjanpalo H: System identifiability based on the power series expansion of the solution. Mathematical Biosciences. 1978, 41: 2133. 10.1016/00255564(78)900639.View Article
 Walter E, Lecourtier Y: Global approaches to identifiability testing for linear and nonlinear state space models. Mathematics and Computers in Simulation. 1982, 24: 472482. 10.1016/03784754(82)906450.View Article
 Vajda S, Godfrey K, Rabitz H: Similarity transformation approach to identifiability analysis of nonlinear compartmental models. Mathematical Biosciences. 1989, 93: 217248. 10.1016/00255564(89)900242.View ArticlePubMed
 Diop S, Fliess M: On nonlinear observability. First European Control Conference, ECC'91, Grenoble. 1991, 152157.
 Fliess M, Glad T: An algebraic approach to linear and nonlinear control. Essays on Control: Perspectives in the Theory and its Applications. Edited by: Treutelman HL, Willeuis JC. 1993, Boston: Birkhauser, 223267.View Article
 Ljung L, Glad T: On global identifiability of arbitrary model parametrizations. Automatica. 1994, 30: 265276. 10.1016/00051098(94)900299.View Article
 Margaria G, Riccomagno E, White LJ: Structural identifiability analysis of some highly structured families of statespace models using differential algebra. Journal of Mathematical Biology. 2004, 49: 433454. 10.1007/s0028500302613.View ArticlePubMed
 Margaria G, Riccomagno E, Chappell MJ, Wynn HP: Differential algebra methods for the study of the structural identifiability of rational function statespace models in the biosciences. Mathematical Biosciences. 2001, 174: 126. 10.1016/S00255564(01)000797.View ArticlePubMed
 Saccomani M, Audoly S, D'Angio L: Parameter identifiability of nonlinear systems: the role of initial conditions. Automatica. 2003, 39: 619632. 10.1016/S00051098(02)003023.View Article
 Yates J, Evans N, Chappell M: Structural identifiability analysis via symmetries of differential equations. Automatica. 2009, 45: 25852591. 10.1016/j.automatica.2009.07.009.View Article
 Walter E, Lecourtier Y, Happel J: On the structural output distinguishability of parametric models, and its relations with structural identifiability. IEEE Transactions on Automatic Control. 1984, AC29: 5657.View Article
 Walter E, Pronzato L: On the identifiability and distinguishability of nonlinear parametric models. Mathematics and Computers in Simulation. 1996, 42: 125134. 10.1016/03784754(95)001239.View Article
 Evans ND, Chappell MJ, Chapman MJ, Godfrey KR: Structural indistinguishability between uncontrolled (autonomous) nonlinear analytic systems. Automatica. 2004, 40: 19471953. 10.1016/j.automatica.2004.06.002.View Article
 Pohjanpalo H, Wahlström B: On the uniqueness of linear compartmental systems. International Journal of Systems Science. 1977, 8: 619632. 10.1080/00207727708942068.View Article
 Yates J, Jones R, Walker M, Cheung S: Structural identifiability and indistinguishability of compartmental models. Expert Opinion on Drug Metabolism and Toxicology. 2009, 5: 295302. 10.1517/17425250902773426.View ArticlePubMed
 Godfrey K, Chapman M, Vajda S: Identifiability and indistinguishability of nonlinear pharmacokinetic models. Journal of Pharmacokinetics and Pharmacodynamics. 1994, 22: 229251.View Article
 Davidescu F, Jorgensen S: Structural parameter identifiability analysis for dynamic reaction networks. Chemical Engineering Science. 2008, 63: 47544762. 10.1016/j.ces.2008.06.009.View Article
 Vajda S, Rabitz H: Identifiability and distinguishability of general reaction systems. Journal of Physical Chemistry. 1994, 98: 52655271. 10.1021/j100071a016.View Article
 Vajda S, Rabitz H: Identifiability and distinguishability of firstorder reaction systems. Journal of Physical Chemistry. 1988, 92: 701707. 10.1021/j100314a024.View Article
 Craciun G, Pantea C: Identifiability of chemical reaction networks. Journal of Mathematical Chemistry. 2008, 44: 244259. 10.1007/s109100079307x.View Article
 Szederkényi G: Comment on "Identifiability of chemical reaction networks" by G. Craciun and C. Pantea. Journal of Mathematical Chemistry. 2009, 45: 11721174. 10.1007/s1091000894998.View Article
 Pohjanpalo H, Wahlström B: Software for solving identification and identifiability problems, e.g. in compartmental systems. Mathematics and Computers in Simulation. 1982, 24: 490493. 10.1016/03784754(82)906486.View Article
 Raksányi A, Lecourtier Y, Walter E, Venot A: Identifiability and distinguishability testing via computer algebra. Mathematical Biosciences. 1985, 77: 245266. 10.1016/00255564(85)901002.View Article
 Bellu G, Saccomani M, Audoly S, D'Angio L: DAISY: A new software tool to test global identifiability of biological and physiological systems. Computer Methods and Programs in Biomedicine. 2007, 88: 5261. 10.1016/j.cmpb.2007.07.002.PubMed CentralView ArticlePubMed
 Saccomani M, Audoly S, Bellu G, D'Angio L: Examples of testing global identifiability of biological and biomedical models with the DAISY software. Computers in Biology and Medicine. 2010, 40: 402407. 10.1016/j.compbiomed.2010.02.004.PubMed CentralView ArticlePubMed
 Liao J, Boscolo R, Yang Y, Tran L, Sabatti C, Roychowdhury V: Network component analysis: reconstruction of regulatory signals in biological systems. Proceedings of the National Academy of Sciences of the United States of America. 2003, 100: 1552215527. 10.1073/pnas.2136632100.PubMed CentralView ArticlePubMed
 Yue H, Brown M, Knowles J, Wang H, Broomhead D, Kell D: Insights into the behaviour of systems biology models from dynamic sensitivity and identifiability analysis: a case study of an NF κB signalling pathway. Molecular Biosystems. 2006, 2: 640649. 10.1039/b609442b.View ArticlePubMed
 Zak D, Gonye G, Schwaber J, Doyle F: Importance of input perturbations and stochastic gene expression in the reverse engineering of genetic regulatory networks: insights from an identifiability analysis of an in silico network. Genome Research. 2003, 13: 23962405. 10.1101/gr.1198103.PubMed CentralView ArticlePubMed
 Chen WW, Niepel M, Sorger PK: Classic and contemporary approaches to modeling biochemical reactions. Genes & Development. 2010, 24: 18611875. 10.1101/gad.1945410.View Article
 Jaqaman K, Danuser G: Linking data to models: data regression. Nature Reviews Molecular Cell Biology. 2006, 7: 813819. 10.1038/nrm2030.View ArticlePubMed
 Banga JR, BalsaCanto E: Parameter estimation and optimal experimental design. Essays in Biochemistry. 2008, 45: 195209. 10.1042/BSE0450195.View ArticlePubMed
 Ashyraliyev M, FomekongNanfack Y, Kaandorp J, Blom J: Systems biology: parameter estimation for biochemical models. FEBS Journal. 2009, 276: 886902. 10.1111/j.17424658.2008.06844.x.View ArticlePubMed
 Soliman S, Heiner M: A unique transformation from ordinary differential equations to reaction networks. PLoS ONE. 2010, 5: e1428410.1371/journal.pone.0014284. (pp. 16)PubMed CentralView ArticlePubMed
 Werhli AV, Husmeier D: Reconstructing gene regulatory networks with Bayesian networks by combining expression data with multiple sources of prior knowledge. Statistical Applications in Genetics and Molecular Biology. 2007, 6: article no. 15
 Mukherjee S, Speed TP: Network inference using informative priors. Proceedings of the National Academy of Sciences of the United States of America. 2008, 105: 1431314318. 10.1073/pnas.0802272105.PubMed CentralView ArticlePubMed
 Szederkényi G, Hangos KM, Péni T: Maximal and minimal realizations of reaction kinetic systems: computation and properties. MATCH Communications in Mathematical and in Computer Chemistry. 2011, 65: 309332.
 Szederkényi G: Computing sparse and dense realizations of reaction kinetic systems. Journal of Mathematical Chemistry. 2009, 47: 551568.View Article
 Szederkényi G, Hangos KM: Finding complex balanced and detailed balanced realizations of chemical reaction networks. Journal of Mathematical Chemistry. 2011, 49: 11631179. 10.1007/s1091001198049.View Article
 Johnston MD, Siegel D: Linear conjugacy of chemical reaction networks. Journal of Mathematical Chemistry. 2011, 49: 12631282. 10.1007/s1091001198174.View Article
 Wang Y, Zhang XS, Chen L: Optimization meets systems biology. BMC Systems Biology. 2010, 4 (Suppl 2): S110.1186/175205094S2S1. (pp. 14)PubMed CentralView ArticlePubMed
 Dantzig GB, Thapa MN: Linear Programming 1: Introduction. 1997, SpringerVerlag
 Floudas C: Nonlinear and MixedInteger Optimization. 1995, Oxford University Press
 Raman R, Grossmann IE: Integration of logic and heuristic knowledge in MINLP optimization for process synthesis. Computers and Chemical Engineering. 1992, 16: 155171. 10.1016/00981354(92)85003Q.View Article
 Feinberg M: Lectures on chemical reaction networks.Notes of lectures given at the Mathematics Research Center, University of Wisconsin. 1979, [http://www.che.eng.ohiostate.edu/%7efeinberg/LecturesOnReactionNetworks/]
 Hárs V, Tóth J: On the inverse problem of reaction kinetics. Qualitative Theory of Differential Equations. Edited by: Farkas M, Hatvani L. 1981, NorthHolland, Amsterdam, 30: 363379.
 Hangos KM, Szederkényi G: Mass action realizations of reaction kinetic system models on various time scales. Journal of Physics: Conference Series (5th International Workshop on MultiRate Processes and Hysteresis). 2011, 268: 012009
 Jokar S, Pfetsch ME: Exact and approximate sparse solutions of underdetermined linear equations. SIAM Journal on Scientific Computing. 2008, 31: 2344. 10.1137/070686676.View Article
 Donoho DL: For most large undetermined systems of linear equations the minimal l1norm solution is also the sparsest solution. Communications on Pure and Applied Mathematics. 2006, 59: 903934.
 Zavlanos MM, Julius AA, Boyd SP, Pappas GJ: Inferring stable genetic networks from steadystate data. Automatica. 2011, 47: 11131122. 10.1016/j.automatica.2011.02.006.View Article
 Szederkényi G, Hangos KM, Tuza Z: Finding Weakly Reversible Realizations of Chemical Reaction Networks Using Optimization. MATCH Communications in Mathematical and in Computer Chemistry. 2012, 67: 193212.
 Kim H, Bond R: Multicore software technologies: a survey. IEEE Signal Processing Magazine. 2009, 26: 8089.View Article
 Mileyko Y, Joh RI, Weitz JS: Smallscale copy number variation and largescale changes in gene expression. Proceedings of the National Academy of Sciences of the United States of America. 2008, 105: 1665916664. 10.1073/pnas.0806239105.PubMed CentralView ArticlePubMed
 Alon U: An Introduction to Systems Biology: Design Principles of Biological Circuits. 2007, Chapman & Hall/CRC Mathematical & Computational Biology
 Ljung L: Perspectives on system identification. Annual Reviews in Control. 2010, 34: 112. 10.1016/j.arcontrol.2009.12.001.View Article
 Conradi C, Flockerzi D, Raisch J, Stelling J: Subnetwork analysis reveals dynamic features of complex (bio)chemical networks. Proceedings of the National Academy of Sciences of the United States of America. 2007, 104: 1917519180. 10.1073/pnas.0705731104.PubMed CentralView ArticlePubMed
 Müller S, Hofbauer J, Endler L, Flamm C, Widder S, Schuster P: A generalized model of the repressilator. Journal of Mathematical Biology. 2006, 53: 905937. 10.1007/s0028500600359.View ArticlePubMed
 Dilao R, Muraro D: A software tool to model genetic regulatory networks. Applications to the modeling of threshold phenomena and of spatial patterning in drosophila. PLoS ONE. 2010, 5: e1074310.1371/journal.pone.0010743. (pp. 110)PubMed CentralView ArticlePubMed
 Bansal M, di Bernardo D: Inference of gene networks from temporal gene expression profiles. IET Systems Biology. 2007, 1: 306312. 10.1049/ietsyb:20060079.View ArticlePubMed
 BangJensen J, Gutin G: Digraphs: Theory, Algorithms and Applications. 2001, Springer
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.