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 [13]. 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 [513], signal transduction [1417] and metabolic networks [1825].
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. [2830] 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 [3335]. 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 [4143]. 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 [4648]. 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 [5760].
The importance of identifiability has been recognized previously in systems biology, too [14,6164]. However, and despite a number of works illustrating ways to test the structural and practical identifiability of models [6567], 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 [7174], 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) [7578]. The analysis of four biological system models can be found in the 'Results and discussion' section, followed by the conclusions.
Methods
The model class considered in this paper is of the following form
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
Following [79] and several other works, we will characterize CRNs with the following three sets.
1.
2.
where α_{ij }are nonnegative integers and are called the stoichiometric coefficients.
3.
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.
Assuming massaction kinetics, the following dynamical description will be used to describe the timeevolution of species concentrations [29,79]:
where x_{i }denotes the concentration of species X_{i}. Let us denote the transpose and the (i, j)th element of an arbitrary matrix W by W^{T }and W_{i,j}, respectively, where i is the row index and j is the column index. The jth column of Y contains the composition of complex C_{j}, i.e. Y_{i,j }= α_{ji}. The structure and parameters of the reaction graph are stored in the column conservation matrix A_{k }(also called the Kirchhoff matrix of the CRN) as follows
Finally, ψ : ℝ^{n }↦ℝ^{m }is a monomialtype vector mapping defined by
Dynamical equivalence of massaction networks
As it is known even from the early literature [28], CRNs with different structures and/or parametrization can give rise to the same
kinetic differential equations. Therefore, we will call two CRNs given by the matrix
pairs
where for
In this case,
We will assume throughout the paper that the set of complexes (i.e. the stoichiometric matrix Y) is fixed and known before the computations. In this case, the condition (6) for dynamical equivalence can be written as
where
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 [7173,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
The following is a straightforward extension of the results published in [71]. To prove the uniqueness of a CRN structure given a set of simple constraints, we
have to extend the notions of dense and sparse realizations. The constraint set denoted
by
where s is the number of individual constraints, and i_{k }≠ j_{k }for k = 1, ..., s. Now we can introduce the following definitions. A dynamically equivalent constrained realization of a CRN (Y, A_{k}) is a reaction network
P1 Given a CRN (Y, A_{k}) and a constraint set
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.
The task of determining which reactions of a CRN belong to the dense realization can be effectively solved through the following problem set consisting of m(m  1) LP computation steps, where m is the number of complexes in the CRN.
where the decision variables are the offdiagonal entries of A_{k}, and U_{ij }are appropriately large positive upper bounds for [A_{k}]_{i,j }to exclude the possibility of unbounded feasible solutions. The reaction C_{q }→ C_{p }is in the dense realization if and only if the maximal objective function value for
f_{pq }in (10) is positive. Let us denote the solution of (10) corresponding to (p, q), p ≠ q by
By construction, ε_{ij }≥ 0 ∀i ≠ j, and ε_{ij }> 0 if and only if the reaction
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
We will call a reaction a core reaction, if it is present in any dynamically equivalent realization of a CRN with a given complex set (and possibly an additional constraint set). Other reactions, the rate coefficient of which can be zero in certain realizations, are called noncore reactions. It clearly follows from the definition, but is remarked separately that the set of core reactions is generally not identical to the set of reactions of a sparse realization. The identification of core reactions of a CRN has not been published yet, therefore we give the outline of the corresponding computation method. Firstly, a dense realization of the network has to be computed to get all the mathematically possible reactions. Then, for each reaction C_{p }→ C_{q }in the dense realization, the feasibility of the following constraint set has to be checked:
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
Let us recall eq. (1). Shortly speaking, global structural identifiability means that
where
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
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.
Additional file 1. Detailed numerical data of the CRNs shown in Examples 13. This file contains the detailed data (i.e. stoichiometric matrices and reaction rate coefficients) of the dynamically equivalent reaction networks studied in Examples 1,2 and 3. The individual sheets correspond to the different examples.
Format: XLS Size: 120KB Download file
This file can be viewed with: Microsoft Excel Viewer
Example 1: a positive feedback motif
The first example is a positive feedback motif shown in Figure 1a and taken from [87] containing 5 species, 11 complexes and 9 reactions. This basic motif is also discussed in [88]. The network contains a gene that promotes its own transcription and translation after dimerization. In the model, X_{1 }and X_{2 }denote the concentrations of protein monomers and dimers, respectively. X_{3 }and X_{4 }are the concentrations of unoccupied and occupied promoters, respectively, and X_{5 }corresponds to the mRNA. The degradation of dimers is ignored. The roles of the reaction rate coefficients are the following: k_{1 }and k_{2 }are the dimerization and redimerization rates, respectively. k_{3 }and k_{4 }are the binding and dissociation rates of the dimer to the promoter, while k_{5 }and k_{6 }denote the activated and basal transcription rates, respectively. k_{7 }is the degradation rate of the mRNA, k_{8 }is the degradation rate of the monomer, and k_{9 }denotes the translation rate. The timeevolution of the speciesconcentrations is described by the following ODEs:
Figure 1. Positive feedback motif: original reaction graph and dense realization structure. (a) This subfigure shows the reaction graph of a gene regulation network model with positive feedback described originally in [87] and used in Example 1. (b) This subfigure shows all the mathematically possible reactions that can result in the same dynamical behaviour as the original biologically meaningful network shown in Figure 1. The corereactions in the dense realization are shown with solid arrows, while the noncore reactions are indicated by dashed arrows.
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.
The matrices characterizing the stoichiometry and graph structure of the system are the following (indicating only the nonzero nondiagonal elements of A_{k}):
We used the following parameter values that were taken from the Appendix of [87].
where the units of measure are [M^{1}] for k_{1}, ..., k_{4}, and [min^{1}] for k_{5}, ..., k_{9}. The dynamically equivalent dense realization of the network is shown in Figure 1b, where the 8 core and 4 noncore reactions are indicated separately. The three different sparse structures are shown in the subplots of Figure 2. The first subplot is identical to the original structure shown in Figure 1a. This means that the mechanism cannot be described exactly with less than 9 reactions. It turns out from the second and third subplots that (at least mathematically), the degradation of mRNA is dynamically not a necessary element of the model. However, the biological plausibility of the mathematically possible structures and reactions always has to be carefully examined.
Figure 2. Sparse realization structures for the positive feedback motif. Three different dynamically equivalent structures can be given for the positive feedback motif with the minimal number of reactions. The core and noncore reactions are indicated in the same way as in Figure 1.b.
As it is expected, the possible structures of sparse/dense realizations and the corresponding core and noncore reactions can change with the modification of parameter values. This is illustrated in Figure 3a, where the following randomly generated parameter values were used:
Figure 3. The effect of modifying the complex set and the parameters. (a) The core and noncore reactions of the dense realization of the positive feedback motif are shown in this subfigure with a randomly selected parametrization that is different from the one given in [87]. (b) The core and noncore reactions of the dense realization of the positive feedback motif can be seen in this subfigure when an additional complex X_{2 }+ X_{4 }is involved into the model.
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.
In the next step, let us assume that another complex, namely X_{2 }+ X_{4 }is allowed in the model (again not necessarily assuming biological meaningfulness in this particular case). With the addition of this new complex, the stoichiometric matrix of the system can be written as
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
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:
x_{1}: [Sic1], x_{2}: [Sic1P], x_{3}: [Clb], x_{4}: [Clb·Sic1], x_{5}: [Clb·Sic1P], x_{6}: [Cdc14], x_{7}: [Sic1P·Cdc14], x_{8}: [Clb·Sic1P·Cdc14], x_{9}: [Clb·Sic1·Clb]. The original structure with 18 reactions is shown in Figure 4a. The Y matrix of the network is given by
Figure 4. Model of a biochemical switch in yeast cells. (a) The subfigure shows the original structure of a CRN describing a biochemical switch published in [90]. The numbering of species and rate coefficients is identical to the description in the original paper. (b) The dense realization of the network is depicted in this subfigure and contains 28 reactions, out of which only 12 belong to the set of core reactions.
The nonzero offdiagonal elements of A_{k }are (the diagonal ones can be computed using the column conservation property):
Since there are no parameter values published in [90], we used the following randomly selected rate coefficients:
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
Consider the repressilator model shown in Figure 5 with 5 nodes where also autoactivation is assumed. Similarly to [91], we make the following assumptions: cooperative regulator binding, genes are present in constant amounts, transcription and translation are modeled by singlestep kinetics, and finally, proteins are degraded by first order reactions. We note that complex dynamic phenomena such as multiple steady states or oscillations have been shown with a wide range of parameters in similar systems, especially in the case when the number of genes is odd [91]. We also assume that there is some protein production (leakage) when both the activator and the repressor are bound to the genes (although this assumption does not affect the main results of the forthcoming analysis). It is clearly shown in [92] that kinetic models with simple massaction kinetics very effectively describe complex dynamics in genetic regulatory networks, therefore we follow the same modeling methodology. Using the assumptions listed above, the CRN describing the system is the following:
Figure 5. A standard repressilator structure. A repressilator structure with 5 nodes and autoactivation is shown in the figure. The massaction type CRN model of this structure contains 51 distinct complexes and 55 reactions.
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).
Two cases with different sets of randomly selected rate coefficients were studied, and the structures of the obtained results were the same. The numerical details can be found on the 3rd sheet of Additional File 1: CRN_data.xls. The total number of reactions for the repressilator model is 55 that is equal to the number of reactions in the sparse realization. The dense realization contains 70 reactions which means that there are a maximum of 15 more mathematically possible reactions while maintaining exactly the same dynamics as the original biological model. These additional reactions are the following:
The number of core reactions in the model are 45. The set of noncore reactions (that, in principle can be substituted by other reactions) is given by
In particular, it is easy to show (see also Additional File 1: CRN_data.xls) that reactions
Example 4: sparse linear gene regulation network models
For structural identification, gene regulation networks are often modeled as linear timeinvariant systems [84,93] of the form
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.
First, consider the 'true' genetic network structure that was simulated and inferred in [93] and that is redrawn in Figure 6a. From the figure, we can reconstruct the structure of the corresponding A matrix as follows (the exact parameter values are not described in the paper, but the investigated structural properties do not depend on the individual parameter values)
Figure 6. A sparse gene regulation network and their structural identifiability properties. (a) This subfigure is the reproduction of one of the sparse gene regulation networks used for structural identification in [93]. The network has 11 activation (solid edges), 6 suppression (dashed edges), and 3 autoregulation links (at nodes 1, 3 and 6) with undefined sign. (b) The subfigure shows that the vast majority of the randomly generated models in the case when N = 2 and K = 2 are structurally unidentifiable, because not all nodes of the network are reachable from the perturbed one. (c) As the network becomes less sparse (N = 10, K = 3), the structural identifiability properties are quickly improving. In this case, more than 80% of the randomly generated models are structurally identifiable.
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.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
All authors contributed to the conception and design of the work. JRB and GS selected and evaluated the examples. GS performed the numerical computations. All authors contributed to the writing of the manuscript. All authors read and approved the final manuscript.
A Appendix
Proof of P1. Let us denote the ith column of any matrix W by W_{·,i}. The proof is based on the following well known fact of linear algebra. Consider an inhomogeneous set of linear equations:
If x = p is any particular solution of (44) then the entire solution set for (44) can be characterized as
The matrix equation Y · A_{k }= M (see eqs. (3) and (8)) obviously defines m sets of linear equations of the form
Let us choose any i indexing the sets of equations in (46). For simplicity, let p = [A_{k}]_{·,i}, b = M_{·,i}. Let us assume that there are z elements of the constraint set (9) where j_{k }= i for k = 1, ..., s. (If z is 0, then we get the earlier result proved in [71].) These constraints can be expressed by further linear equations of the form:
The equation sets (46) and (47) can be written into a single set of equations as
where
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.
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.
References

Wolkenhauer O: Systems biology: The reincarnation of systems theory applied in biology?
Briefings in Bioinformatics 2001, 2:258270. PubMed Abstract  Publisher Full Text

Stelling J: Mathematical models in microbial systems biology.
Current Opinion in Microbiology 2004, 7:513518. PubMed Abstract  Publisher Full Text

Kitano H: Computational systems biology.
Nature 2002, 420:206210. PubMed Abstract  Publisher Full Text

Csete M, Doyle J: Reverse engineering of biological complexity.
Science 2002, 295:16641669. PubMed Abstract  Publisher Full Text

De Jong H: Modeling and simulation of genetic regulatory systems: a literature review.
Journal of Computational Biology 2002, 9:67103. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

de la Fuente A, Brazhnik P, Mendes P: Linking the genes: inferring quantitative gene networks from microarray data.
Trends in Genetics 2002, 18:395398. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Schaber J, Klipp E: Modelbased inference of biochemical parameters and dynamic properties of microbial signal transduction networks.
Current Opinion in Biotechnology 2010, 22:109116. PubMed Abstract  Publisher Full Text

Li S, Assmann S, Albert R: Predicting essential components of signal transduction networks: a dynamic model of guard cell abscisic acid signaling.

SaezRodriguez J, Kremling A, Conzelmann H, Bettenbrock K, Gilles E: Modular analysis of signal transduction networks.

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:56.
(pp. 126)
PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text 
Arkin A, Shen P, Ross J: A test case of correlation metric construction of a reaction pathway from measurements.
Science 1997, 277:12751279. Publisher Full Text

Stelling J, Klamt S, Bettenbrock K, Schuster S, Gilles E: Metabolic network structure determines key aspects of functionality and regulation.
Nature 2002, 420:190193. PubMed Abstract  Publisher Full Text

Jeong H, Tombor B, Albert R, Oltvai Z, Barabási A: The largescale organization of metabolic networks.
Nature 2000, 407:651654. PubMed Abstract  Publisher Full Text

Forster J, Famili I, Fu P, Palsson B, Nielsen J: Genomescale reconstruction of the Saccharomyces cerevisiae metabolic network.
Genome Research 2003, 13:244. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kell D: Metabolomics and systems biology: making sense of the soup.
Current Opinion in Microbiology 2004, 7:296307. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text

Yuan Y, Stan GB, Warnick S, Goncalves J: Robust dynamical network structure reconstruction.
Automatica 2011, 47:12301235. Publisher Full Text

Stolovitzky G, Monroe D, Califano A: Dialogue on ReverseEngineering Assessment and Methods.
Annals of the New York Academy of Sciences 2007, 1115:122. PubMed Abstract  Publisher Full Text

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:e9202.
(pp. 118)
PubMed Abstract  Publisher Full Text  PubMed Central Full Text 
Horn F, Jackson R: General mass action kinetics.
Archive for Rational Mechanics and Analysis 1972, 47:81116.

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. Publisher Full Text

Érdi P, Tóth J: Mathematical Models of Chemical Reactions. In Theory and Applications of Deterministic and Stochastic Models. Manchester, Princeton: Manchester University Press, Princeton University Press; 1989.

Ljung L: System Identification: Theory for the User. 2nd edition. Prentice Hall; 1999.

Walter E, Pronzato L: Identification of Parametric Models. Springer; 1997.

Bellmann R, Aström KJ: On structural identifiability.
Mathematical Biosciences 1970, 7:329339. Publisher Full Text

Walter E: Identification of State Space Models. Springer; 1982.

Walter E: Identifiability of Parametric models. Pergamon Press, Oxford; 1987.

Pohjanpalo H: System identifiability based on the power series expansion of the solution.
Mathematical Biosciences 1978, 41:2133. Publisher Full Text

Walter E, Lecourtier Y: Global approaches to identifiability testing for linear and nonlinear state space models.
Mathematics and Computers in Simulation 1982, 24:472482. Publisher Full Text

Vajda S, Godfrey K, Rabitz H: Similarity transformation approach to identifiability analysis of nonlinear compartmental models.
Mathematical Biosciences 1989, 93:217248. PubMed Abstract  Publisher Full Text

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. In Essays on Control: Perspectives in the Theory and its Applications. Edited by Treutelman HL, Willeuis JC. Boston: Birkhauser; 1993:223267.

Ljung L, Glad T: On global identifiability of arbitrary model parametrizations.
Automatica 1994, 30:265276. Publisher Full Text

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. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text

Saccomani M, Audoly S, D'Angio L: Parameter identifiability of nonlinear systems: the role of initial conditions.
Automatica 2003, 39:619632. Publisher Full Text

Yates J, Evans N, Chappell M: Structural identifiability analysis via symmetries of differential equations.
Automatica 2009, 45:25852591. Publisher Full Text

Walter E, Lecourtier Y, Happel J: On the structural output distinguishability of parametric models, and its relations with structural identifiability.

Walter E, Pronzato L: On the identifiability and distinguishability of nonlinear parametric models.
Mathematics and Computers in Simulation 1996, 42:125134. Publisher Full Text

Evans ND, Chappell MJ, Chapman MJ, Godfrey KR: Structural indistinguishability between uncontrolled (autonomous) nonlinear analytic systems.
Automatica 2004, 40:19471953. Publisher Full Text

Pohjanpalo H, Wahlström B: On the uniqueness of linear compartmental systems.
International Journal of Systems Science 1977, 8:619632. Publisher Full Text

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. PubMed Abstract  Publisher Full Text

Godfrey K, Chapman M, Vajda S: Identifiability and indistinguishability of nonlinear pharmacokinetic models.
Journal of Pharmacokinetics and Pharmacodynamics 1994, 22:229251.

Davidescu F, Jorgensen S: Structural parameter identifiability analysis for dynamic reaction networks.
Chemical Engineering Science 2008, 63:47544762. Publisher Full Text

Vajda S, Rabitz H: Identifiability and distinguishability of general reaction systems.
Journal of Physical Chemistry 1994, 98:52655271. Publisher Full Text

Vajda S, Rabitz H: Identifiability and distinguishability of firstorder reaction systems.
Journal of Physical Chemistry 1988, 92:701707. Publisher Full Text

Craciun G, Pantea C: Identifiability of chemical reaction networks.
Journal of Mathematical Chemistry 2008, 44:244259. Publisher Full Text

Szederkényi G: Comment on "Identifiability of chemical reaction networks" by G. Craciun and C. Pantea.
Journal of Mathematical Chemistry 2009, 45:11721174. Publisher Full Text

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. Publisher Full Text

Raksányi A, Lecourtier Y, Walter E, Venot A: Identifiability and distinguishability testing via computer algebra.
Mathematical Biosciences 1985, 77:245266. Publisher Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chen WW, Niepel M, Sorger PK: Classic and contemporary approaches to modeling biochemical reactions.
Genes & Development 2010, 24:18611875. PubMed Abstract  Publisher Full Text

Jaqaman K, Danuser G: Linking data to models: data regression.
Nature Reviews Molecular Cell Biology 2006, 7:813819. PubMed Abstract  Publisher Full Text

Banga JR, BalsaCanto E: Parameter estimation and optimal experimental design.
Essays in Biochemistry 2008, 45:195209. PubMed Abstract  Publisher Full Text

Ashyraliyev M, FomekongNanfack Y, Kaandorp J, Blom J: Systems biology: parameter estimation for biochemical models.
FEBS Journal 2009, 276:886902. PubMed Abstract  Publisher Full Text

Soliman S, Heiner M: A unique transformation from ordinary differential equations to reaction networks.
PLoS ONE 2010, 5:e14284.
(pp. 16)
PubMed Abstract  Publisher Full Text  PubMed Central Full Text 
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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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.

Szederkényi G, Hangos KM: Finding complex balanced and detailed balanced realizations of chemical reaction networks.
Journal of Mathematical Chemistry 2011, 49:11631179. Publisher Full Text

Johnston MD, Siegel D: Linear conjugacy of chemical reaction networks.
Journal of Mathematical Chemistry 2011, 49:12631282. Publisher Full Text

Wang Y, Zhang XS, Chen L: Optimization meets systems biology.
BMC Systems Biology 2010, 4(Suppl 2):S1.
(pp. 14)
PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text 
Dantzig GB, Thapa MN: Linear Programming 1: Introduction. SpringerVerlag; 1997.

Floudas C: Nonlinear and MixedInteger Optimization. Oxford University Press; 1995.

Raman R, Grossmann IE: Integration of logic and heuristic knowledge in MINLP optimization for process synthesis.
Computers and Chemical Engineering 1992, 16:155171. Publisher Full Text

Feinberg M: [http://www.che.eng.ohiostate.edu/%7efeinberg/LecturesOnReactionNetworks/] webcite
Lectures on chemical reaction networks.Notes of lectures given at the Mathematics Research Center, University of Wisconsin. 1979.

Hárs V, Tóth J: On the inverse problem of reaction kinetics. In Qualitative Theory of Differential Equations. Volume 30. Edited by Farkas M, Hatvani L. NorthHolland, Amsterdam; 1981::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. Publisher Full Text

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. Publisher Full Text

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.

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Alon U: An Introduction to Systems Biology: Design Principles of Biological Circuits. Chapman & Hall/CRC Mathematical & Computational Biology; 2007.

Ljung L: Perspectives on system identification.
Annual Reviews in Control 2010, 34:112. Publisher Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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. PubMed Abstract  Publisher Full Text

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:e10743.
(pp. 110)
PubMed Abstract  Publisher Full Text  PubMed Central Full Text 
Bansal M, di Bernardo D: Inference of gene networks from temporal gene expression profiles.
IET Systems Biology 2007, 1:306312. PubMed Abstract  Publisher Full Text

BangJensen J, Gutin G: Digraphs: Theory, Algorithms and Applications. Springer; 2001.