Email updates

Keep up to date with the latest news and content from BMC Systems Biology and BioMed Central.

Open Access Methodology article

Analyzing fixed points of intracellular regulation networks with interrelated feedback topology

Nicole Radde

Author Affiliations

Institute for Systems Theory and Automatic Control University of Stuttgart Pfaffenwaldring 9, Germany

BMC Systems Biology 2012, 6:57  doi:10.1186/1752-0509-6-57

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1752-0509/6/57


Received:12 December 2011
Accepted:6 June 2012
Published:6 June 2012

© 2012 Radde; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Modeling the dynamics of intracellular regulation networks by systems of ordinary differential equations has become a standard method in systems biology, and it has been shown that the behavior of these networks is often tightly connected to the network topology. We have recently introduced the circuit-breaking algorithm, a method that uses the network topology to construct a one-dimensional circuit-characteristic of the system. It was shown that this characteristic can be used for an efficient calculation of the system’s fixed points.

Results

Here we extend previous work and show several connections between the circuit-characteristic and the stability of fixed points. In particular, we derive a sufficient condition on the characteristic for a fixed point to be unstable for certain graph structures and demonstrate that the characteristic does not contain the information to decide whether a fixed point is asymptotically stable. All statements are illustrated on biological network models.

Conclusions

Single feedback circuits and their role for complex dynamic behavior of biological networks have extensively been investigated, but a transfer of most of these concepts to more complex topologies is difficult. In this context, our algorithm is a powerful new approach for the analysis of regulation networks that goes beyond single isolated feedback circuits.

Keywords:
Circuit-breaking algorithm; feedback circuit; fixed point analysis; fixed point bifurcation

Background

Describing the dynamic behavior of molecular interactions in a cell or cell compartment by chemical reaction kinetics has become a standard approach in systems biology for metabolic pathways as well as for regulatory networks. Since qualitative knowledge about these interactions is often available from experiments, literature or databases, which can be represented as network graphs, several different graph-based approaches have been developed to analyze the behavior of the networks. These methods operate solely on the graphs without detailed knowledge of the kinetic rates. They show for example that certain subnetwork structures are necessary to generate complex behavior such as oscillations, hysteresis or multistationarity. Thus, such behavior can be excluded for relatively small and simple networks that lack these subnetworks. So far, most of these approaches have the following limitations for practical use: First, they only allow to make statements for relatively simple graph topologies, and second, they are often restricted to very specific model classes such as metabolic networks of the form <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M1">View MathML</a> with stoichiometric matrix S and (often polynomial) flux vector v(x) [1] or regulatory networks whose Jacobian matrices have constant signs on the off-diagonal elements [2-5]. Similar analysis methods that work for more complex graph topologies or more general network model classes are rare. On the other hand, it has been shown in various contexts that interrelated feedback structures contribute to the robustness of intracellular regulation processes [6-13]. In most studies this is demonstrated by analyzing a specific model via simulations with varying parameter values, for example via Monte Carlo simulations. Although the conclusions from these studies are very helpful and valuable, it is not clear to which extend they can be generalized to other network models. These results, which show the importance of feedback in regulation processes, provide a further basis for the need of new methods that can deal with interrelated feedback in dynamic network models in a more general way. We expect that the more complex the graph topologies, the more does the system’s behavior depend on the kinetic rate laws, and less can be concluded from the structure alone. Thus, these new methods can probably not be completely independent of equations and parameters.

A new approach in this direction has been introduced in our previous work [14] for a general class of regulatory network models. We introduced the circuit-breaking algorithm (CBA), a method which operates on the graph topology to construct a one-dimensional characteristic that inherits important information about the behavior of the system. In particular, we demonstrated that the zeros of this characteristic are related to the system’s fixed points.

In this paper we extend this work and show that the characteristic contains information about stability of the fixed points and can furthermore be used to detect bifurcation point candidates. The paper is structured as follows: We give a brief overview over our network model class and the circuit-breaking algorithm and show how it works on a network for cellular differentiation of hematopoietic stem cells [15]. Based on these results, we investigate relations between the stability of fixed points and the slope of the circuit-characteristic that is constructed by the CBA. It is shown that a negative slope at a zero of the characteristic does generally not contain any information about the stability of the respective fixed point, while a positive slope implies that the fixed point is unstable, at least for certain graph topologies. We demonstrate results on biological network models for tryptophan regulation in Escherichia coli[11] and the repressilator model [16].

Results and Discussion

The circuit-breaking algorithm

Here we introduce the regulatory network model class and summarize the concept of the CBA. For details we refer to [14]. Since the formal description of the algorithm is very technical and needs a lot of indices, we will thereafter directly show how it works on a concrete network example, from which we hope that it makes the general concept easier understandable.

We consider regulatory networks models that are described by a system of first order ordinary differential equations

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M2">View MathML</a>

(1)

with underlying interaction graph (I-graph) G(V ,E) that illustrates the dependence structure of the variables, i.e.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M3">View MathML</a>

(2)

and

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M4">View MathML</a>

(3)

Trajectories of the system should be bounded, a biologically plausible assumption which also implies that the system has at least one fixed point. It is sometimes useful for the analysis to introduce sign-labels of edges in the I-graph if the corresponding partial derivative is monotone, which means that the influence of a component on another one is purely activating or purely inhibiting regardless of the state of the system. Contrary to many other methods, the CBA does not rely on this monotonicity assumption.

Given a regulatory network model, i.e. a differential equation system <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M5">View MathML</a> and the I-graph topology G(V ,E), the CBA consists of the following steps:

1. Find strongly connected components of G(V,E):The first step of the CBA is a partitioning of the vertex set V into strongly connected components (SCC), i.e. the maximal strongly connected subgraphs, which we denote by Gk(VkEk), k=1,…,K. The new graph GSCC(VSCCESCC), which is obtained by shrinking all vertices of a SCC to one single vertex and drawing an edge <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M6">View MathML</a> between two vertices <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M7">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M8">View MathML</a> of this graph when there exist vertices <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M9">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M10">View MathML</a> that are connected with an edge <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M11">View MathML</a> in the original graph G(VE), has a hierarchical topology without any circuits. This fact is illustrated in Figure 1.We numerate the SCCs according to this hierarchical order in the network. Fixed point coordinates of the system can iteratively be calculated for each SCC, starting with the SCC on top and following the hierarchical structure downwards. In this scheme the fixed point coordinates of the SCCs upstream serve as constant input u for subsequent SCCs. An example for this concept of iterative fixed point calculation for SCCs is given in [14]. We denote these sets of fixed points of SCC k with input u by <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M12">View MathML</a>. For the sake of simplicity we skip the index u in the following, but bear in mind that the fixed point set Fkhas to be calculated for each input u.

2. Construct characteristics for each SCC in dependence of input u and calculate the fixed points from it’s zeros: The core of the CBA is the construction of a one-dimensional characteristic <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M13">View MathML</a> for a SCC Gk(Vk,Ek) for each input u. This is done in the following way:

(a) Find the set C of all elementary circuits and list them as set of vertex subsets

(b) Find a minimal circuit-covering vertex set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M14">View MathML</a> such that at least one element of each subset in C is contained in <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M15">View MathML</a> and set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M16">View MathML</a>. Collect the rest of the vertices in the set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M17">View MathML</a>. Relabel vertices such that <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M18">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M19">View MathML</a>.

(c) Break all circuits by removing all edges that point to vertices of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M20">View MathML</a>. Mathematically, this is done by setting these variables to fixed input values κ=(κ1,…,κm), i.e. xi=κi. The result is an acyclic or hierarchical graph topology.

(d) The fixed point coordinates of variables in <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M21">View MathML</a>, denoted by <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M22">View MathML</a>, can be calculated in dependence of these inputs κ.

(e) The circuits are iteratively closed by releasing the vertices in <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M23">View MathML</a> one after another, starting backwards with vm. This translates into shifting the respective vertex vi from the set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M24">View MathML</a> to <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M25">View MathML</a>, reducing the vector κby the same element, and solving the implicit equation of the form

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M26">View MathML</a>

(4)

for xi to get the set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M27">View MathML</a> of fixed point coordinates of the variable xi in dependence of the vector κ. The set F(κ) has to be updated accordingly. Equation (4) has to be solved numerically. For i=2,…,mwe denote the expression on the left hand side of equation (4) partial circuit-characteristic. The number of input variables of these characteristics is reduced by one in each step, since κis reduced by one element in each step. Thus, when releasing the last vertex v1 in <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M28">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M29">View MathML</a> is a one-dimensional characteristic that is called the circuit-characteristicc(κ1) of the respective SCC. It’s zeros correspond to the fixed point coordinates of variable x1, denoted by <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M30">View MathML</a>. If we leave the current SCC and go to the next one, we refer to this characteristic as <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M31">View MathML</a>.

(f) The corresponding fixed point coordinates of the other variables can be calculated iteratively by inserting the values of the set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M32">View MathML</a> into the partial circuit-characteristics in reverse order. These coordinates are then collected in the set F of fixed point coordinates of the SCC k. If we leave the SCC k, we refer to this set as Fk.

thumbnailFigure 1. Division of vertices into strongly connected components. An example of a graph G(V ,E) and it’s partitioning into strongly connected components (left). Within a strongly connected component, each pair of vertices is connected via a path. If each SCC is contracted to a single vertex, the resulting graph is circuit-free (right) (if it contained circuits, the SCCs are not maximal, since all SCCs in the circuit could be merged to a larger SCC) and thus has a hierarchical order.

The structure of the CBA is illustrated in Figure 2 with a flow chart.

thumbnailFigure 2. Flow chart of the Circuit-Breaking Algorithm. Flow chart of the Circuit-Breaking Algorithm. Blue boxes indicate that these calculations are done within a SCC of the graph, yellow boxes describe the iterative closing of the circuits within this SCC by releasing vertices in the set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M33">View MathML</a> one after another. The green boxes refer to actions on the full graph G(V ,E).

Application of the CBA to a model for hematopoietic stem cell differentiation

To motivate the subsequent investigations on the characteristics of regulatory network models and it’s relation to fixed points and their stability, we consider a network model for the cellular differentiation of hematopoietic stem cells described in [15]:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M34">View MathML</a>

(5)

This model describes the interplay between the two lineage-specific counter-acting suppressors Gfi-1 (x2) and Egr(1,2) (x3) during cellular differentiation for the neutrophil and macrophage cell fate choices, respectively. These are activated by their transcription factors C/EBPα(x1) and PU.1 (x4), respectively. Together, they regulate the expression of lineage-specific downstream genes, which are not further specified in the model and denoted by Mac (x5) and Neut (x6). The model was build based on chemical reaction kinetics that describe interaction of the molecular species. The cellular state is assumed to be directly correlated to the fixed point concentrations of the transcription factors, as described further below. Furthermore, the model was non-dimensionalized after some simplifications by rescaling time and protein concentrations. The two parameters that are left, eN and eM, are the rescaled synthesis rate of the transcription factor C/EBPα, which is not regulated in the model, and the maximal rescaled synthesis rate of the transcription factor PU.1.

Figure 3 shows the bifurcation diagram of all six variables with bifurcation parameter μ=eM and condition eN=eM that was created using xpppaut. For eM=0 the system has a globally stable fixed point at x=0. The system undergoes a saddle-node bifurcation at <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M35">View MathML</a>. It has a globally stable fixed point for <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M36">View MathML</a> and two stable fixed points divided by an intermediate unstable one for <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M37">View MathML</a>. It can also be seen that the stable fixed point branch that exists for all eMrepresents the neutrophil state, since the fixed point coordinates of the neutrophil specific proteins (x1,x2,x6) increase monotonically along this branch. The macrophage state is represented by the stable fixed point branch that appears at <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M38">View MathML</a>.

thumbnailFigure 3. Bifurcation diagrams of a model for cellular differentiation of hematopoietic stem cells. Bifurcation diagrams of the hematopoietic stem cell differentiation network (system (5)) described in [15] with bifurcation parameter μ=eM. Has been created with Additional file 1.

Now we use the CBA to construct the characteristic of this system and compare this with the information of the bifurcation diagram. As can be seen in Figure 4, the I-graph of system (5) consists of four strongly connected components given by <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M39">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M40">View MathML</a> with circuit sets <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M41">View MathML</a>, and minimal circuit-covering vertex sets <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M42">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M43">View MathML</a>.

thumbnailFigure 4. I-graph of a model for cellular differentiation of hematopoietic stem cells. I-graph G(V ,E) of system (5). It consists of four SCCs, as indicated with the grey boxes, and C2has two interrelated positive circuits.

We start with G1(V1,E1), which does not contain any circuits. Thus, we just have to solve <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M44">View MathML</a> in system (5), which leads to the set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M45">View MathML</a> of fixed points of G1. The fixed point set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M46">View MathML</a> of G2(V2,E2) is calculated by breaking the two circuits at x2, i.e. setting <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M47">View MathML</a>. Inserting <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M48">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M49">View MathML</a> into <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M50">View MathML</a> leads to the circuit-characteristic

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M51">View MathML</a>

(6)

which can, by inserting the respective terms for the synthesis rates r, be rewritten as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M52">View MathML</a>

(7)

with

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M53">View MathML</a>

(8)

This characteristic is shown in Figure 5 (center row), along with the sets <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M54">View MathML</a>, i=3,4,5,6, for parameter values eM={0.2,0.3221,0.5} (left, center, right row, respectively).

The following properties of the system can be identified from these figures:

1. The fixed point coordinates of all variables x3,x4,x5and x6behave monotonically with the input κ2, which represents the neutrophil state. The macrophage specific proteins x3,x4and x5decrease with increasing <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M55">View MathML</a>, x6increases.

2. Looking at the characteristics (center row) for different values eM, it is monotonically decreasing for <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M56">View MathML</a> (left), and thus has a single zero, which corresponds to the single fixed point branch for <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M57">View MathML</a>. For the value eM=0.2, which is chosen here, we get the fixed point <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M58">View MathML</a>={0.2,0.81,0.43,0.14,0.86,1.75}, as indicated in the graphs. This state represents an intermediate non-differentiated progenitor cell state. The saddle-node bifurcation is represented by the second zero of the characteristic that appears at <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M59">View MathML</a> (center column). The respective fixed point set is <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M60">View MathML</a>0.59} and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M61">View MathML</a>0.10,0.17,2.22}.Finally, the characteristic has three zeros for <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M62">View MathML</a> (right column) and thus the system has three fixed points in this range. For the chosen value eM=0.5 we can read off the fixed point set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M63">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M64">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M65">View MathML</a>. Here, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M66">View MathML</a> represents the macrophage state, where Egr and PU.1 are highly expressed, and C/EBPαis low, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M67">View MathML</a> stands for the neutrophil state in which C/EBPαis dominant, and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M68">View MathML</a> is an unstable intermediate state that separates the two basins of attraction.

thumbnailFigure 5. Circuit characteristics of a model for cellular differentiation of hematopoietic stem cells. From top to bottom: Sets <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M69">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M70">View MathML</a>, circuit-characteristic <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M71">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M72">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M73">View MathML</a> for values eM=0.2 (left column), <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M74">View MathML</a> (center column) and eM=0.5 (right column) of system (5). The fixed points that correspond to the zeros of the characteristic are also indicated. Has been created by Additional file 2.

Seeking for further parallels between the bifurcation diagram (Figure 3) and the characteristics in Figure 5, the question arises if the characteristic also contains information about bifurcations and stability of the fixed points. Clearly, the parameters for which the characteristic touches the x-axis without intersection are bifurcation value candidates. Furthermore, looking at this example, a self-evident guess would be to assume that stability can be determined in the same way as for one-dimensional vector fields: The fixed points are stable if the slope of the characteristic at the corresponding zero κ is negative, i.e. <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M75">View MathML</a>, and it is unstable if the slope is positive, i.e. <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M76">View MathML</a>. We will further investigate these assumptions in the following subsections. In order to do so, we consider in the following strongly connected I-graphs, which allows to neglect the indices u and k, such that indexing can be simplified. The results are, however, easily transferable to arbitrary graphs, since construction of the characteristic is done separately for each strongly connected component. We will continue by denoting the characteristic simply with c(κ1), where <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M77">View MathML</a> is the value of the variable x1, the one which is released lastly. We first prove the following proposition, which relates the slope of the characteristic to the determinant of the Jacobian matrix Jf(x) of the system:

Proposition 1

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M78">View MathML</a>

(9)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M79">View MathML</a>

(10)

with F(x1) denoting the fixed point coordinates of variables x2,…,x|V|in dependence of x1, and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M80">View MathML</a> is the Jacobian matrix of the subnetwork spanned by the vertices V∖{v1}.

The proof is given in the Methods section. Note that Proposition 1 holds for all inputs κ1, but we are here especially interested in the zeros of the characteristic, i.e. the set of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M81">View MathML</a> with <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M82">View MathML</a>, and we will in the following subsection sometimes denote the corresponding fixed point with <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M83">View MathML</a>, if appropriate.

Instability of fixed points

From Proposition 1 it follows that a positive slope <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M84">View MathML</a> implies that <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M85">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M86">View MathML</a> have the same signs. According to the Hartman-Grobman theorem (see e. g. [17]), a fixed point <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M87">View MathML</a> is unstable if <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M88">View MathML</a> has at least one eigenvalue with positive real part. Unfortunately, we are not aware of a relation between the determinant of Jf(x) and it’s minors that can be used to show the following: If <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M89">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M90','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M90">View MathML</a> have the same signs, this implies the existence of an eigenvalue with positive real part and hence implies instability of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M91">View MathML</a>. Thus we will concentrate on certain graph structures which we call leading vertex graphs (LVG). LVGs are strongly connected components with minimal circuit covering vertex set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M92','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M92">View MathML</a> that consists of one single element v1. In other words, G(V E) has a vertex that is contained in all elementary circuits, and hence the characteristic c(κ1) can be constructed in a single circuit-closing step. The I-graph of the hematopoietic stem cell differentiation network consists of SCCs that are all LVGs, while the two networks considered in the proof of proposition (1) do not belong to this class, because two circuit-closing steps were necessary in each of these cases. For LVGs we can show that a positive slope of the characteristic at a zero implies instability of the respective fixed point. The proof is given in the Methods section.

Stability of fixed points

In contrast to a positive slope, a negative slope of the circuit-characteristic at a fixed point coordinate κ alone does not contain information about the stability of the respective fixed point. We demonstrate this with two examples. The first one consists of a single negative feedback circuit, the repressilator model described in [16]. This is a synthetic transcriptional network of the three repressor proteins lacI, tetR and cI and their corresponding promoters, which was constructed to create periodic expression in Escherichia coli:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M93">View MathML</a>

(11)

with i={lacI,tetR,cI}, j={cI,lacI,tetR}, and mi and pi are mRNA and protein concentrations, respectively. The system has a trapping region, that is, a positively invariant region in the state space that is eventually reached by all trajectories, which guarantees the existence of at least one fixed point. Bounds are given by <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M94','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M94">View MathML</a><a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M95','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M95">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M96">View MathML</a>i=1,2,3. The I-graph (Figure 6) is strongly connected, the circuit set C consists of one subset that contains all six nodes, C={{mipi}i=1,2,3}, and hence the set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M97">View MathML</a> has one single element and the graph is a LVG.

thumbnailFigure 6. I-graph of the repressilator model. I-graph of the repressilator model (11) described in [16].

Note that because of the symmetry of the model, the circuit-characteristic is independent of the choice of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M98','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M98">View MathML</a> here. It is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M99">View MathML</a>

(12)

which can be simplified to

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M100','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M100">View MathML</a>

(13)

where we have used <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M101','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M101">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M102','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M102">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M103','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M103">View MathML</a>. Equation (13) is strictly decreasing, and, importantly, independent of the parameter β.

On the contrary, the eigenvalues of the Jacobian matrix of the system and hence the stability of the fixed point are not (see also the stability diagram in Figure 1b in [16]). This dependence is illustrated in Figure 7, where the real and imaginary parts of the eigenvalues λ(β) of the Jacobian matrix <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M104','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M104">View MathML</a> are plotted as functions of β for parameter values α=290, n=2 and α0=10. For these parameter values the system has a fixed point <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M105','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M105">View MathML</a> for all i=1,2,3 (that is independent of β). It can be seen that there exist solutions with positive real part for small values of β, and hence the fixed point is unstable in this range. It becomes stable through a Hopf bifurcation for increasing values of β. Thus we have shown that <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M106','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M106">View MathML</a> and in particular the stability of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M107','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M107">View MathML</a> depend on β, while c(κ1) does not. From this example we conclude that our assumption is not true for zeros of the characteristic with negative slope. The corresponding fixed point of the system can generally be stable or unstable. In the Methods section proposition 1 is verified for this example.

thumbnailFigure 7. Eigenvalues of the repressilator model. Real (green) and imaginary (red) parts of the eigenvalues λ(β) of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M108','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M108">View MathML</a> of the repressilator model (11) with parameters α=290, n=2, α0=10 and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M109','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M109">View MathML</a>. The figure was created by calculating the characteristic polynomial χ(λ,β) of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M110','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M110">View MathML</a>, which is given here as <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M111','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M111">View MathML</a> with <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M112">View MathML</a>, and using a Newton gradient search with several random starting points to find the eigenvalues λwith accuracy <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M113','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M113">View MathML</a><10−4. Has been created with Additional files 3 and 4.

As a further example we consider the tryptophan regulation network in Escherichia coli described in [11], which can be written as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M114','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M114">View MathML</a>

(14)

where the state vector x corresponds to the free operator sites (OR), mRNA (M), enzyme (E) and tryptophan (T) concentrations. C(xKm) are sigmoidally decreasing functions,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M115','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M115">View MathML</a>

(15)

This model describes the regulation of the tryptophan concentration via different mechanisms, i.e. genetic regulation via binding of tryptophan to it’s operator site, described by C(x4t1m1), mRNA attenuation (C(x4t2m2)) and enzyme inhibition (C(x4t3m3)). The parameters k1k2k3and k4represent kinetic rate constants for synthesis of free operator, mRNA transcription, translation and tryptophan synthesis, respcetively, K are half-saturation constants for the inhibition processes, Ot denotes the total operator site concentration, and γand μrefer to degradation and diluation rates due to cell growth. The term <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M116','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M116">View MathML</a> describes the uptake of tryptophan for protein synthesis in the cell.

Analyzing this system with the parameter values given in [11] (<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M117">View MathML</a><a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M118','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M118">View MathML</a><a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M119','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M119">View MathML</a><a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M120','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M120">View MathML</a>g=25μMKg=0.2μM) using xppaut and choosing the dilution rate μas bifurcation parameter, the system shows a Hopf bubble between <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M121','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M121">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M122','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M122">View MathML</a> (Figure 8). The system has a unique fixed point that is unstable between these two values and shows sustained oscillations in this range. Outside the Hopf bubble the oscillations are damped and the fixed point is globally stable.

thumbnailFigure 8. Bifurcation diagrams of a model for tryptophan regulation inEscherichia coli. Bifurcation diagrams of the tryptophan regulation model (14) in Escherichia coli described in [11] with dilution rate μas bifurcation parameter. The system shows two Hopf bifurcations at <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M123','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M123">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M124">View MathML</a>. Has been created with Additional file 5.

The corresponding I-graph is shown in Figure 9. It is strongly connected.

The circuit set C and the minimal circuit covering vertex set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M125','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M125">View MathML</a> are <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M126','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M126">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M127">View MathML</a>. Since <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M128">View MathML</a> consists of a single element, this system is a LVG and only one circuit-closing step is necessary to calculate the set of fixed points. The circuit-characteristic can be calculated analytically here and is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M129','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M129">View MathML</a>

(16)

where r4can iteratively be calculated via

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M130">View MathML</a>

(17)

As can be seen in Figure 10, c(κ4) is strictly decreasing (bottom row), since all circuits in the graph are negative.

thumbnailFigure 9. I-graph of the tryptophan regulation network. I-graph of the tryptophan regulation network (14).

Furthermore, the fixed points of the system can be determined by the zeros of the characteristic, as depicted in the figure: For μ=0.01, for example, c(κ4) has a zero at <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M131','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M131">View MathML</a>, which corresponds to the fixed point coordinate <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M132','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M132">View MathML</a>. Inserting this value into <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M133','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M133">View MathML</a>, i=1,2,3, we get the fixed point <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M134','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M134">View MathML</a>, and likewise for the other dilution rates. The qualitative courses of c(κ4) and also for the fixed point sets <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M135','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M135">View MathML</a> do not differ for the three dilution rates. In particular, the slope of the characteristic is in all three cases negative at the zero. However, the bifurcation diagrams in Figure 3 indicate that the respective fixed points are stable for μ=0.01 and μ=0.2, but unstable for μ=0.1. Thus this is a further example that a negative slope of the characteristic at a zero does not imply stability of the respective fixed point.

thumbnailFigure 10. Circuit-characteristics of the tryptophan regulation model. Sets <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M136','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M136">View MathML</a> for i=1,2,3 and circuit-characteristic c(κ4) of the model of tryptophan regulation in Escherichia coli with the same parameter set as was used in [11] and dilution parameters μ=0.01 (right column), μ=0.1 (center column) and μ=0.2 (left column). Has been created with Additional file 6.

Conclusions

In this paper we have extended previous work on the analysis of fixed points for regulatory network models. Based on the circuit-breaking algorithm, which was introduced in [14] and which uses the topology of the interaction graph to construct a one-dimensional circuit-characteristic whose zeros correspond to the fixed points of the system, we further investigated this characteristic with respect to fixed points of the system and their stability. Here we demonstrated that the characteristic is in some aspects similar to a one-dimensional vector field and that the CBA is also useful to find fixed point bifurcations. Information about the stability of fixed points can partly be derived from the slopes at the respective zeros of the characteristic. We used our methods to analyze the fixed points of models for hematopoietic stem cell differentiation, tryptophan regulation in Escherichia coli and the repressilator in Escherichia coli. In particular, we have shown that a positive slope of the characteristic at a zero can imply instability, at least for certain graph topologies, which we call leading vertex graphs. These are characterized by leading vertices for all strongly connected components that are contained in all circuits. Although we have noticed that many network models belong to this model class, this restriction on the topology for sure limits the use of our approach. However, we believe that the implication can further be generalized to other network topologies, although a pure translation of the techniques that we are currently using is not possible. Thus a generalization is one topic for future work.

On the contrary, generally no conclusions about stability can be drawn from a negative slope, and the respective fixed point can either be stable or unstable. If it is unstable, we interpret this result as a kind of time-delay. This delay is due to the response time of the network to changes in the input κi. It is not visible in the characteristic any more, where the effects of all feedback circuits have been summarized to a single effective one comprising only one component. This effect might be similar to a time-delay that destabilizes a stable fixed point in a one-dimensional vector field.

While this manuscript was in revision, we became aware of a recent paper [18] that seems to be closely related to our work in some aspects. In this paper, small phosphorylation motifs in signaling pathways are investigated subject to their ability to show bistable behavior. The authors follow the same idea of variable elimination to construct finally one-dimensional functions that contain information about the fixed points of the system and their stability. However, the techniques used therein are build on mass action kinetics and rational functions and explicitly use mass conservation relations. However, some of the mathematical ideas behind that seem to be related to our work, and a further comparison would be interesting.

Generally, the efficiency of the CBA and the analysis introduced here depends on the graph topology and the complexity to solve the implicit equations therein. Construction of the circuit-characteristic is particularly simple and efficient for graph topologies whose strongly connected components have minimal circuit-covering vertex set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M137','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M137">View MathML</a> with only few elements, and thus our theory can be particularly helpful to analyze such networks.

In the future we will try to generalize results further, such that our approach is applicable to a broad range of regulatory network models. We will also further investigate the connection between the partial circuit-characteristics and the influence of the respective sets of circuits that are closed on the coordinates, number and stability of the system’s fixed points. We believe that our analysis can lead to the identification of circuit sets which are responsible for certain behaviors of the system that are connected to bifurcations of fixed points. Finally, we hope that we can contribute towards developing analysis methods that facilitate an understanding of the role of interrelated feedback circuits in regulatory network models for the system’s overall behavior.

Methods

In this section we collect the mathematical technicalities that are needed to show the statements made in the Results and Discussion section of the manuscript.

Proof of Proposition 1

This section shows the proof of Proposition 1. To avoid complex indexing, the relation is exemplarily shown on a fully connected 3-vertex network and a network with four vertices. These examples are non-trivial in the sense that the cardinality of the minimal covering vertex set, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M138','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M138">View MathML</a>, contains more than one element, such that calculation of the characteristic requires more than a single circuit-closing step. Thus the principles of these two examples can be generalized to other I-graphs.

3-vertex model

Proof

We consider a regulatory network model with a fully-connected I-graph with three vertices:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M139','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M139">View MathML</a>

(18)

whose Jacobian matrix is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M140','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M140">View MathML</a>

(19)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M141','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M141">View MathML</a>

(20)

We now construct the circuit-characteristic c(κ1) using the CBA, whose steps are illustrated in Figure 11.

thumbnailFigure 11. Circuit-breaking algorithm for a regulatory network model with three vertices. The circuit-breaking algorithm for a regulatory network model with three vertices and fully connected I-graph.

In order to calculate it’s derivative and show Proposition 1, we will repeatedly use the Implicit function theorem (IFT), which reads:

Implicit function theorem[19]: Let U be an open set in <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M142','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M142">View MathML</a> and let <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M143','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M143">View MathML</a> be a <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M144','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M144">View MathML</a> function with k≥1. Consider a point <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M145','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M145">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M146','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M146">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M147','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M147">View MathML</a>, with <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M148','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M148">View MathML</a>. If the n×n matrix <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M149','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M149">View MathML</a> of partial derivatives is invertible, then there are open sets <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M150','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M150">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M151','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M151">View MathML</a> with <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M152','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M152">View MathML</a> and a unique Ck function ψ:VmVn such that f(xψ(x))=c for all xVm. Moreover, f(xy)≠c if (xy)∈Vm×Vn and yψ(x). The derivative of the function ψis given by the formula

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M153','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M153">View MathML</a>

(21)

In the first step we break all circuits by fixing x1=κ1and x2=κ2(Figure 11left) and get the partial circuit-characteristic and the fixed point set

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M154','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M154">View MathML</a>

(22)

with derivative given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M155','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M155">View MathML</a>

(23)

Here we have used the IFT with m=2, n=1, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M156','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M156">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M157','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M157">View MathML</a>, c=0, and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M158','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M158">View MathML</a>.

In the next step we release v2(Figure 11center) and get the partial circuit-characteristic and the fixed point set

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M159','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M159">View MathML</a>

(24)

with derivative

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M160','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M160">View MathML</a>

(25)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M161','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M161">View MathML</a>

(26)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M162','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M162">View MathML</a>

(27)

Here we have used the IFT with m=1, n=1, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M163','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M163">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M164','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M164">View MathML</a>, c=0, and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M165','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M165">View MathML</a>.

In the last step also vertex v1 is released (Figure 11right). The circuit-characteristic c(x1) reads:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M166','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M166">View MathML</a>

(28)

and its derivative is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M167','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M167">View MathML</a>

(29)

Multiplying this expression with <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M168','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M168">View MathML</a>leads to

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M169','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M169">View MathML</a>

(30)

4-vertex model

Proof

Additionally, we outline the proof of proposition 1 for a non-trivial four-component network (Figure 12):

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M170','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M170">View MathML</a>

(31)

thumbnailFigure 12. Circuit-breaking algorithm for a regulatory network model with four vertices. The circuit-breaking algorithm for a regulatory network model with four vertices.

whose Jacobian matrix is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M171','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M171">View MathML</a>

(32)

where we used the same notation as before, i.e. <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M172','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M172">View MathML</a>.

Again we construct the circuit-characteristic c(κ1) using the CBA and the IFT for it’s derivatives. First we break all circuits by fixing x1=κ1 and x2=κ2(Figure 12left) and calculating the fixed point coordinates of the remaining vertices:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M173','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M173">View MathML</a>

(33)

with derivatives

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M174','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M174">View MathML</a>

(34)

and

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M175','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M175">View MathML</a>

(35)

with derivatives

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M176','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M176">View MathML</a>

(36)

and

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M177','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M177">View MathML</a>

(37)

In the next step we release v2(Figure 12center) and get the partial circuit-characteristic and the fixed point set

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M178','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M178">View MathML</a>

(38)

whose derivative is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M179','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M179">View MathML</a>

(39)

Thus we have expressed the fixed point coordinates of x3 and x4 in terms of κ1, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M180','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M180">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M181','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M181">View MathML</a>. Finally we release v1. The circuit-characteristic c1(x1) reads:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M182','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M182">View MathML</a>

(40)

and its derivative is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M183','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M183">View MathML</a>

(41)

Using again the IFT to eliminate derivatives of fixed point coordinates, i.e.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M184','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M184">View MathML</a>

(42)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M185','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M185">View MathML</a>

(43)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M186','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M186">View MathML</a>

(44)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M187','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M187">View MathML</a>

(45)

the derivative of the characteristic becomes

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M188','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M188">View MathML</a>

(46)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M189','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M189">View MathML</a>

(47)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M190','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M190">View MathML</a>

(48)

Setting <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M191','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M191">View MathML</a>, we can see that equation (49) equals <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M192','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M192">View MathML</a> when multiplied with βf33f44. This can easily be seen by multiplying (49) out and rearranging the order of the summands. □

Unstable fixed points in LVGs

Since for LVGs the subnetwork spanned by the set V∖{v1} does by definition not contain any circuits, we get a simple expression for <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M193','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M193">View MathML</a>, namely:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M194','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M194">View MathML</a>

(49)

with <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M195','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M195">View MathML</a> for all i. Thus, the sign of this expression is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M196','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M196">View MathML</a>

(50)

Now we assume that <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M197','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M197">View MathML</a> is stable, and hence R(λ)<0 for all eigenvalues λof <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M198','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M198">View MathML</a>. It follows that

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M199','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M199">View MathML</a>

(51)

In any case, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M200','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M200">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M201','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M201">View MathML</a> have different signs, a contradiction to <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M202','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M202">View MathML</a>, which completes the proof.

Verification of Proposition 1 for the repressilator model

Let us verify Proposition 1 for the repressilator model. We can identify

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M203','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M203">View MathML</a>

(52)

and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M204','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M204">View MathML</a>, such that both <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M205','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M205">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M206','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M206">View MathML</a> depend on β, but c(κ1) does not, since β is canceled out. The derivative of c(κ1) (equation (13)) with respect to κ1reads

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M207','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M207">View MathML</a>

(53)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M208','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M208">View MathML</a>

(54)

which equals <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M209','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M209">View MathML</a>.

Endnotes

1This plot was created with the program xppaut and Additional file 7 in the Supplement (tryptophanmodel.ode)

2This plot was generated with the program gnuplot and Additional file 8 (tryptophanmodel-c.gp) in the Supplement.

Competing interests

The authors declare that they have no competing interests.

Author’s contributions

NR designed research and carried out calculations on the examples.

Additional files

Addtional file 1. Stemcellmodel. This file was used to create the bifurcation diagrams of the hematopoietic stem cell model using the program xppaut (Figure 3).

Format: ODE Size: 1KB Download fileOpen Data

Addtional file 2. Stemcellmodel-02. The file stemcellmodel-02.gp was used to create the circuit-characteristic c(κ) and fixed point sets F(κ) of the hematopoietic stem cell model (Figure 5) with bifurcation parameter μ=0.2 using the program gnuplot.

Format: GP Size: 3KB Download fileOpen Data

Addtional file 3. Stemcellmodel-bif. The file stemcellmodel-bif.gp was used to create the circuit-characteristic c(κ) and fixed point sets F(κ) of the hematopoietic stem cell model (Figure 5) with bifurcation parameter μ=μusing the program gnuplot.

Format: GP Size: 4KB Download fileOpen Data

Addtional file 4. Stemcellmodel-05. The file stemcellmodel-05.gp was used to create the circuit-characteristic c(κ) and fixed point sets F(κ) of the hematopoietic stem cell model (Figure 5) with bifurcation parameter μ=0.5 using the program gnuplot.

Format: GP Size: 4KB Download fileOpen Data

Addtional file 5. Newton-2d. This python script was used to calculate the eigenvalues of the Jacobian matrix <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M210','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M210">View MathML</a> of the repressilator model at it’s fixed point using a Newton gradient search with random starting points. These eigenvalues were written into the file ‘eigenvalues.txt’.

Format: PY Size: 5KB Download fileOpen Data

Addtional file 6. Eigenvalues. This file contains the eigenvalues of the Jacobian matrix <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/57/mathml/M211','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/57/mathml/M211">View MathML</a> of the repressilator model at it’s fixed point and was created by running the program newton-2d.py. It was used to create Figure 7.

Format: TXT Size: 196KB Download fileOpen Data

Addtional file 7. Tryptophanmodel. This file was used to create the bifurcation diagrams of the tryptophan regulation model using the program xppaut (Figure 8).

Format: ODE Size: 1KB Download fileOpen Data

Addtional file 8. Tryptophanmodel-c. This file was used to create the circuit-characteristic c(κ) and fixed point sets F(κ) of the tryptophan regulation model (Figure 10) using the program gnuplot.

Format: GP Size: 4KB Download fileOpen Data

Acknowledgements

This work was supported by the German Research Foundation (DFG) within the Cluster of Excellence in Simulation Technology (EXC 310/1) and the funding programme Open Access Publishing.

References

  1. Conradi C, Flockerzi D, Raisch J, Stelling J: Subnetwork analysis reveals dynamic features of complex (bio)chemical networks.

    Proc Natl Acad Sci USA 2007, 104(49):19175. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  2. Gouzé JL: Positive and negative circuits in dynamical systems.

    J Biol Syst 1998, 6(21):11. OpenURL

  3. Radde N, Bar N, Banaji M: Graphical methods for analysing feedback in biological networks - A survey.

    Int J Syst Sci 2010, 41:35. Publisher Full Text OpenURL

  4. Thieffry D: Dynamical roles of biological regulatory circuits.

    Brief Bioinform 2007, 8(4):220. PubMed Abstract | Publisher Full Text OpenURL

  5. Thomas R: On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations. In Numerical methods in the study of critical phenomena, Volume 9 of Springer series. Edited by Della-Dora J, Demongeot J, Lacolle B. Springer series; 1981:180-193. OpenURL

  6. Brunner M, Káldi K: Interlocked feedback loops of the circadian clock of Neurospora crassa.

    Mol Microbiol 2008, 68(2):255. PubMed Abstract | Publisher Full Text OpenURL

  7. Clodong S, Dühring U, Kronk L, Wilde A, Axmann I, Herzel H, Kollmann M: Functioning and robustness of a bacterial circadian clock.

    Mol Syst Biol 2007, 3(90):1. OpenURL

  8. Kwon YK, Cho KH: Quantitative analysis of robustness and fragility in biological networks based on feedback dynamics.

    Bioinformatics 2008, 24(7):987. PubMed Abstract | Publisher Full Text OpenURL

  9. Ma W, Lai L, Ouyang Q, Tang C: Robustness and modular design of the Drosophila segment polarity network.

    Mol Syst Biol 2006., 2(70) OpenURL

  10. Radde N: The role of feedback mechanisms in biological network models - A tutorial.

    Asian J Control 2011, 13(5):597. Publisher Full Text OpenURL

  11. Venkatesh K, Bhartiya S, Ruhela A: Multiple feedback loops are key to a robust dynamic performance of tryptophan regulation in Escherichia coli.

    FEBS Letters 2004, 563:234. PubMed Abstract | Publisher Full Text OpenURL

  12. Wagner A: Circuit topology and the evolution of robustness in two-gene circadian oscillators.

    Proc Natl Acad Sci USA 2005, 102(33):11775. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  13. Wang R, Chen L, Aihara K: Detection of cellular rhythms and global stability within interlocked feedback systems.

    Math Biosci 2007, 209:171. PubMed Abstract | Publisher Full Text OpenURL

  14. Radde N: Fixed point characterization of differential equations with complex graph topology.

    Bioinformatics 2010, 26(22):2874. PubMed Abstract | Publisher Full Text OpenURL

  15. Laslo P, et al.: Multilineage transcriptional priming and determination of alternate hematopoietic cell fates.

    Cell 2006, 126:755. PubMed Abstract | Publisher Full Text OpenURL

  16. Elowitz M, Leibler S: A synthetic oscillatory network of transcriptional regulators.

    Nature 2000, 403:335. PubMed Abstract | Publisher Full Text OpenURL

  17. Guckenheimer J, Holmes P: Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. No. 42 in Applied Mathematical Sciences, Springer series, Springer Verlag, New York, Berlin Heidelberg, Tokyo; 1990. OpenURL

  18. Feliu E, Wiuf C: Enzyme sharing as a cause of multistationarity in signaling systems.

    J R Soc Interface 7 2012, 9(71):1224. Publisher Full Text OpenURL

  19. Hale J, Kocak H: Dynamics and Bifurcations. Texts in Applied Mathematics; OpenURL