Inference of complex biological networks: distinguishability issues and optimization-based solutions

Abstract

Background

The inference of biological networks from high-throughput 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 high-throughput 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 network-inference (sometimes called reverse-engineering [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 one-size-fits-all 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 gene-network 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 noise-free measurements and error-free modeling [33–35]. One of the early approaches for identifiability testing of nonlinear models is based on the Taylor-series 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 Volterra-series coefficients of the system which is the nonlinear generalization of the Laplace-transform 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 Lie-algebras [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 mass-action 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 time-evolution of species concentrations [28, 55]. A related problem, namely the non-unique 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 optimization-based 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 so-called 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 non-uniqueness 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

The model class considered in this paper is of the following form

$$\begin{array}{c}ẋ=f\left(x,u,\theta \right),x\left(\mathsf{\text{0}}\right)=x_0\\ y=h\left(x,u,\theta \right),\end{array}$$

(1)

where x ∈ ℝ ⁿ 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 mass-action 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 mass-action networks

Following [79] and several other works, we will characterize CRNs with the following three sets.

• 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,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 'Feinberg-Horn-Jackson 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₁, C₂, ... 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₁, ..., 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 mass-action kinetics, the following dynamical description will be used to describe the time-evolution of species concentrations [29, 79]:

$$ẋ=Y\cdot A_k\cdot \psi \left(x\right),$$

(3)

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

$${\left[A_k\right]}_{i,j}=\left\{\begin{array}{cc}\hfill -{\sum }_{l=1,l\ne i}^m{k}_{il},\hfill & \hfill \mathsf{\text{if}}i=j\hfill \\ \hfill k_{ji},\hfill & \hfill \mathsf{\text{if}}i\ne j.\hfill \end{array}\right.$$

(4)

Finally, ψ : ℝ ⁿ ↦ℝ ^m is a monomial-type vector mapping defined by

$${\psi }_j\left(x\right)=\prod _{i=1}^n{x}_i^{Y_{i,j}},j=1,\dots ,m.$$

(5)

Dynamical equivalence of mass-action 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 $\left(Y^{\left(1\right)},A_k^{\left(1\right)}\right)$ and $\left(Y^{\left(2\right)},A_k^{\left(2\right)}\right)$dynamically equivalent, if

$$Y^{\left(1\right)}{A}_k^{\left(1\right)}{\psi }^{\left(1\right)}\left(x\right)=Y^{\left(2\right)}{A}_k^{\left(2\right)}{\psi }^{\left(2\right)}\left(x\right)=f\left(x\right),\forall x\in {\overline{ℝ}}_{+}^n,$$

(6)

where for $i=1,2,Y^{\left(i\right)}\in {ℝ}^{n\times m_i}$ have nonnegative integer entries, $A_k^{\left(i\right)}$ are valid Kirchhoff matrices, and

$${\psi }_j^{\left(i\right)}\left(x\right)=\prod _{k=1}^n{x}_k^{{\left[Y^{\left(i\right)}\right]}_{k,j}},i=1,2,j=1,\dots ,m_i.$$

(7)

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.

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

$$Y\cdot A_k^{\left(1\right)}=Y\cdot A_k^{\left(2\right)}=:M,$$

(8)

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 non-unique (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 super-structure) [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 non-diagonal 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 NP-hard 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 MILP-based 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 $\mathcal{K}$ will be used for the exclusion of selected reactions from the CRN, i.e. it is of the form:

$$\mathcal{K}=\left\{{\left[A_k\right]}_{i_1,j_1}=0,\dots ,{\left[A_k\right]}_{i_s,j_s}=0\right\},$$

(9)

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 $\left(Y,A_k^{\prime }\right)$ such that $Y\cdot A_k=Y\cdot A_k^{\prime }$ and the prescribed constraints $\mathcal{K}$ in the form of eq. (9) are fulfilled for $A_k^{\prime }$. A dynamically equivalent constrained dense realization of a CRN (Y, A _k ) is a constrained realization that contains the maximal number of nonzero elements in $A_k^{\prime }$. Similarly, the constrained sparse realization is a constrained realization with the minimal number of nonzeros in $A_k^{\prime }$. To characterize constrained dense/sparse realizations, the results of [71] can be adapted easily as follows.

• 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 MILP-framework in [72]. However, using the structural uniqueness of such realizations given by P1, it is easy to give a polynomial-time 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 super-structure. This allows us to formulate a polynomial-time 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.

$$\begin{array}{c}\mathsf{\text{for each }}p,q=1,\dots ,m,p\ne q\mathsf{\text{do}}:\\ \mathsf{\text{maximize }}{f}_{pq}={\left[A_k\right]}_{p,q}\\ \mathsf{\text{subject to}}:\\ Y\cdot A_k=M,\\ \sum _{i=1}^m{\left[A_k\right]}_{i,j}=0,j=1,\dots ,m,\\ 0\le {\left[A_k\right]}_{i,j}\le U_{ij},i,j=1,\dots ,m,i\ne j,\\ {\left[A_k\right]}_{i,i}\le 0,i=1,\dots ,m,\end{array}$$

(10)

where the decision variables are the off-diagonal 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 ${\overline{A}}_k^{pq}$. Since the linear equality and inequality constraints in (10) are trivially convex, we will use the average of the obtained solutions ${\overline{A}}_k^{pq}$ as a lower bound to compute a possible dense realization in the final optimization step. For this, we define

$${\epsilon }_{ij}={\left[\begin{array}{cc}\hfill \frac{1}{m\left(m-1\right)}\hfill & \hfill \sum _{\begin{array}{c}p,q=1\\ p\ne q\end{array}}^m{Ā}_k^{pq}\hfill \end{array}\right]}_{i,j},i\ne j.$$

(11)

By construction, ε _ij ≥ 0 ∀i ≠ j, and ε _ij > 0 if and only if the reaction ${\mathcal{C}}_j\to {\mathcal{C}}_i$ is in the dense realization. Then the actual dense realization can be determined by solving the following LP feasibility problem for A _k (with arbitrary linear objective function):

$$\begin{array}{c}Y\cdot A_k=M,\\ \sum _{i=1}^m{\left[A_k\right]}_{i,j}=0,j=1,\dots ,m,\\ {\epsilon }_{ij}\le {\left[A_k\right]}_{i,j}\le U_{ij},i,j=1,\dots ,m,i\ne j,\\ {\left[A_k\right]}_{i,i}\le 0,i=1,\dots ,m.\end{array}$$

(12)

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

$$Y\cdot A_k=M$$

(13)

$$\sum _{i=1}^m{\left[A_k\right]}_{i,j}=0,j=1,\dots ,m$$

(14)

$${\left[A_k\right]}_{i,j}\ge 0,i,j=1,\dots ,m,i\ne j,\left(i,j\right)\ne \left(q,p\right)$$

(15)

$${\left[A_k\right]}_{i,i}\le 0,i=1,\dots ,m$$

(16)

$${\left[A_k\right]}_{q,p}=0,$$

(17)

where the matrix A _k contains the decision variables, and the known matrices are Y and M. It is well-known 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 core-reaction 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 multi-core hardware environment [86].

Basic concepts on structural identifiability and distinguishability

Let us recall eq. (1). Shortly speaking, global structural identifiability means that

$$ŷ\left(t|{\theta }^{\prime }\right)\equiv ŷ\left(t|{\theta }^{″}\right)\Rightarrow {\theta }^{\prime }={\theta }^{″},$$

(18)

where

$$ŷ\left(t|\theta \right)=h\left(x\left(t,\theta \right),u\left(t\right),\theta \right),$$

(19)

and x(t, θ) denotes the solution of (1) with parameter vector θ. According to (18), a structurally non-identifiable 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 θ₁ (possibly except for a finite number of values) there is no θ₂ such that the input-output 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 input-output 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 off-diagonal 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 non-identifiability (since zero and non-zero 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 1-3 are contained in a standard spreadsheet form with brief explanations in Additional file 1: CRN_data.xls.

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₁ and X₂ denote the concentrations of protein monomers and dimers, respectively. X₃ and X₄ are the concentrations of unoccupied and occupied promoters, respectively, and X₅ corresponds to the mRNA. The degradation of dimers is ignored. The roles of the reaction rate coefficients are the following: k₁ and k₂ are the dimerization and re-dimerization rates, respectively. k₃ and k₄ are the binding and dissociation rates of the dimer to the promoter, while k₅ and k₆ denote the activated and basal transcription rates, respectively. k₇ is the degradation rate of the mRNA, k₈ is the degradation rate of the monomer, and k₉ denotes the translation rate. The time-evolution of the species-concentrations 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 core-reactions in the dense realization are shown with solid arrows, while the non-core reactions are indicated by dashed arrows.

$${ẋ}_1=-2k_1{x}_1^2+2k_2{x}_2+k_9{x}_5-k_8{x}_1$$

(20)

$${ẋ}_2=k_1{x}_1^2-k_2{x}_2-k_3{x}_2{x}_3+k_4{x}_4$$

(21)

$${ẋ}_3=-k_3{x}_2{x}_3+k_4{x}_4$$

(22)

$${ẋ}_4=k_3{x}_2{x}_3-k_4{x}_4$$

(23)

$${ẋ}_5=k_5{x}_4+k_6{x}_3-k_7{x}_5.$$

(24)

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 well-defined 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 non-diagonal elements of A _k ):

$$Y=\left[\begin{array}{ccccc}2& 0& 0& 0& 0\\ 1& 0& 0& 0& 1\\ 1& 0& 0& 0& 0\\ 0& 1& 1& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 1\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 1\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0\end{array}\right]$$

(25)

$$\begin{array}{c}{\left[A_k\right]}_{2,1}=k_1,{\left[A_k\right]}_{1,2}=k_2,{\left[A_k\right]}_{4,3}=k_3,\\ {\left[A_k\right]}_{3,4}=k_4,{\left[A_k\right]}_{5,4}=k_5,{\left[A_k\right]}_{7,6}=k_6,\\ {\left[A_k\right]}_{9,8}=k_7,{\left[A_k\right]}_{9,10}=k_8,{\left[A_k\right]}_{11,8}=k_9.\end{array}$$

(26)

We used the following parameter values that were taken from the Appendix of [87].

$$\begin{array}{c}{k}_1=k_2=k_3=k_4=10^7,k_5=1.7,\\ k_6=0.025,k_7=0.1,k_8=0.05,k_9=0.5,\end{array}$$

(27)

where the units of measure are [M^-1] for k₁, ..., k₄, and [min^-1] for k₅, ..., k₉. The dynamically equivalent dense realization of the network is shown in Figure 1b, where the 8 core and 4 non-core 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 non-core 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 non-core 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 non-core 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 non-core reactions of the dense realization of the positive feedback motif can be seen in this subfigure when an additional complex X₂ + X₄ is involved into the model.

$$\begin{array}{c}{k}_1=18.9,k_2=7.1,k_3=15.4,\\ k_4=12.7,k_5=10.6,k_6=3.5,\\ k_7=11.3,k_8=9.1,k_9=4.0.\end{array}$$

(28)

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 non-core 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₂ + X₄ 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

$$Y^{\prime }={\left[\begin{array}{ccccc}\hfill 2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \end{array}\right]}^T.$$

(29)

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 non-core 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 non-core 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 17-9\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 mass-action 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₁: [Sic1], x₂: [Sic1P], x₃: [Clb], x₄: [Clb·Sic1], x₅: [Clb·Sic1P], x₆: [Cdc14], x₇: [Sic1P·Cdc14], x₈: [Clb·Sic1P·Cdc14], x₉: [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.

$$Y={\left[\begin{array}{ccccccccc}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T.$$

(30)

The non-zero off-diagonal elements of A _k are (the diagonal ones can be computed using the column conservation property):

$$\begin{array}{c}{\left[A_k\right]}_{2,1}=k_3,{\left[A_k\right]}_{2,3}=k_2,{\left[A_k\right]}_{3,2}=k_1,\\ {\left[A_k\right]}_{4,5}=k_5,{\left[A_k\right]}_{5,4}=k_4,{\left[A_k\right]}_{6,5}=k_6,\\ {\left[A_k\right]}_{6,8}=k_9,{\left[A_k\right]}_{7,8}=k_8,{\left[A_k\right]}_{8,7}=k_7,\\ {\left[A_k\right]}_{9,10}=k_{11},{\left[A_k\right]}_{10,9}=k_{10},{\left[A_k\right]}_{11,10}=k_{12},\\ {\left[A_k\right]}_{12,13}=k_{14},{\left[A_k\right]}_{13,12}=k_{13},{\left[A_k\right]}_{14,13}=k_{15},\\ {\left[A_k\right]}_{15,16}=k_{17},{\left[A_k\right]}_{16,15}=k_{16},{\left[A_k\right]}_{17,16}=k_{18}.\end{array}$$

(31)

Since there are no parameter values published in [90], we used the following randomly selected rate coefficients:

$$\begin{array}{c}k=[4.13.26.77.33.82.44.55.16.2\\ 7.78.69.52.44.95.810.26.38.5{]}^T.\end{array}$$

(32)

The structure of the dense realization indicating the 12 core and 16 non-core 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₅ → X₃ + X₅, X₄ → X₃ + X₄, X₂ + X₃ → X₃ + X₅, and X₃ + X₁ → X₃ + X₄ 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 well-selected reactions leads to the removal of an additional 6 reactions leaving only 18.

Example 3: a repressilator structure with 5 nodes and auto-activation

Consider the repressilator model shown in Figure 5 with 5 nodes where also auto-activation 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 single-step 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 mass-action 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 auto-activation is shown in the figure. The mass-action type CRN model of this structure contains 51 distinct complexes and 55 reactions.

$$G_i+P_i\underset{k_{i,2}}{\overset{k_{i,1}}{\rightleftarrows }}{G}_i^A\left(\mathsf{\text{auto-activation}}1\right)$$

(33)

$$G_i^A\stackrel{k_{i,3}}{\to }{G}_i^A+P_i\left(\mathsf{\text{protein}}\mathsf{\text{production}}1\right)$$

(34)

$$G_i+P_j\underset{k_{i,5}}{\overset{k_{i,4}}{\rightleftarrows }}{G}_i^R\left(\mathsf{\text{repression}}1\right)$$

(35)

$$G_i^R+P_i\underset{k_{i,7}}{\overset{k_{i,6}}{\rightleftarrows }}{G}_i^{AR}\left(\mathsf{\text{auto-activation}}2\right)$$

(36)

$$G_i^A+P_j\underset{k_{i,9}}{\overset{k_{i,8}}{\rightleftarrows }}{G}_i^{AR}\left(\mathsf{\text{repression}}2\right)$$

(37)

$$G_i^{AR}\stackrel{k_{i,10}}{\to }{G}_i^{AR}+P_i\left(\mathsf{\text{protein production}}2\right)$$

(38)

$$P_i\stackrel{k_{i,11}}{\to }0\left(\mathsf{\text{protein degradation}}\right)$$

(39)

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:

$$\begin{array}{c}{G}_i^{AR}\to G_i^R,P_i+G_i^R\to G_i^R,\\ P_i+G_i^R\to P_i+G_i^{AR},\mathsf{\text{for}}i=1,\dots ,5.\end{array}$$

(40)

The number of core reactions in the model are 45. The set of non-core reactions (that, in principle can be substituted by other reactions) is given by

$$G_i^{AR}\leftrightarrows G_i^R+P_i,i=1,\dots ,5.$$

(41)

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

For structural identification, gene regulation networks are often modeled as linear time-invariant systems [84, 93] of the form

$$ẋ=Ax+Bu,$$

(42)

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 auto-activation or auto-repression depending on their sign. x ∈ ℝ ⁿ 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.

$$\left[\begin{array}{cccccccccc}\hfill *\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill +\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill +\hfill & \hfill +\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill *\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill +\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill +\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill +\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill *\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill +\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill +\hfill & \hfill +\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\hfill & \hfill 0\hfill & \hfill -\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\hfill & \hfill +\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill +\hfill & \hfill -\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right],$$

(43)

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 non-identifiability 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 breadth-first-search (DFS or BFS) of the corresponding directed graphs that are fast polynomial-time 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 (3-10) 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 non-identifiable. 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 non-identifiable. 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 non-identifiable. 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 time-complexity 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 polynomial-time 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 MILP-based 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

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:

$$Ax=b$$

(44)

If x = p is any particular solution of (44) then the entire solution set for (44) can be characterized as

$$\mathsf{\text{{}}p+v|v\mathsf{\text{is}}\mathsf{\text{any}}\mathsf{\text{solution}}\mathsf{\text{of}}Ax=\mathsf{\text{0}}}$$

(45)

The matrix equation Y · A _k = M (see eqs. (3) and (8)) obviously defines m sets of linear equations of the form

$$Y\cdot {\left[A_k\right]}_{\cdot ,i}=M_{\cdot ,i},i=1,\dots ,m$$

(46)

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:

$${\left[A_k\right]}_{h,i}=0,h=1,\dots ,z$$

(47)

The equation sets (46) and (47) can be written into a single set of equations as

$$Ȳ\cdot p=\bar{b}$$

(48)

where $Ȳ\in {ℝ}^{\left(n+z\right)\times m}$ and $\bar{b}\in {ℝ}^{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 $Ȳ\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^{″}\ne 0$ for j = 1, ..., l, and there is at least one index l + 1 ≤ c ≤ m for which $p_c^{″}\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.