The circuitbreaking 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
\stackrel{\u0307}{x}=f\left(x\right),\phantom{\rule{1em}{0ex}}x\in {\mathbb{R}}^{n},f\in {\mathcal{C}}^{1}
(1)
with underlying interaction graph (Igraph) G(V ,E) that illustrates the dependence structure of the variables, i.e.
{e}_{j\to i}\in E\phantom{\rule{0.3em}{0ex}}\iff \phantom{\rule{0.3em}{0ex}}\exists \phantom{\rule{0.3em}{0ex}}x\in {\mathbb{R}}^{n}\phantom{\rule{0.3em}{0ex}}\text{such that}\phantom{\rule{0.3em}{0ex}}\frac{\partial {\stackrel{\u0307}{x}}_{i}\left(x\right)}{\partial {x}_{j}}\ne 0
(2)
and
{e}_{i\to i}\in E\phantom{\rule{0.3em}{0ex}}\iff \phantom{\rule{0.3em}{0ex}}\exists \phantom{\rule{0.3em}{0ex}}x\in {\mathbb{R}}^{n}\phantom{\rule{0.3em}{0ex}}\text{such that}\phantom{\rule{0.3em}{0ex}}\frac{\partial {\stackrel{\u0307}{x}}_{i}\left(x\right)}{\partial {x}_{i}}>0.
(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 signlabels of edges in the Igraph 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 \stackrel{\u0307}{x}=f\left(x\right) and the Igraph 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 G^{k}(V^{k}E^{k}), k=1,…,K. The new graph G^{SCC}(V^{SCC}E^{SCC}), which is obtained by shrinking all vertices of a SCC to one single vertex and drawing an edge {e}_{i\to j}^{\mathrm{SCC}} between two vertices {v}_{i}^{\mathrm{SCC}} and {v}_{j}^{\mathrm{SCC}} of this graph when there exist vertices {v}_{i}\in {V}^{i} and {v}_{j}\in {V}^{j} that are connected with an edge {e}_{i\to j} in the original graph G(V E), 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 {}_{u}{F}^{k}. For the sake of simplicity we skip the index u in the following, but bear in mind that the fixed point set F^{k}has 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 onedimensional characteristic {c}^{k}\left({\kappa}_{1}^{k}\right) for a SCC G^{k}(V^{k},E^{k}) 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 circuitcovering vertex set \stackrel{~}{V} such that at least one element of each subset in C is contained in \stackrel{~}{V} and set m=\left\stackrel{~}{V}\right. Collect the rest of the vertices in the set \hat{V}. Relabel vertices such that \stackrel{~}{V}=\{{v}_{1},\dots ,{v}_{m}\} and \hat{V}=\{{v}_{m+1},\dots ,{v}_{\left{V}^{k}\right}\}.

(c)
Break all circuits by removing all edges that point to vertices of \stackrel{~}{V}. Mathematically, this is done by setting these variables to fixed input values κ=(κ_{1},…,κ_{
m
}), i.e. x_{
i
}=κ_{
i
}. The result is an acyclic or hierarchical graph topology.

(d)
The fixed point coordinates of variables in \hat{V}, denoted by F\left(\kappa \right)={\left\{{\stackrel{\u0304}{x}}_{p}\right(\kappa \left)\right\}}_{p=m+1,\dots ,\left{V}^{k}\right}, can be calculated in dependence of these inputs κ.

(e)
The circuits are iteratively closed by releasing the vertices in \stackrel{~}{V} one after another, starting backwards with v_{
m
}. This translates into shifting the respective vertex v_{
i
} from the set \stackrel{~}{V} to \hat{V}, reducing the vector κby the same element, and solving the implicit equation of the form
{f}_{i}({x}_{i},\kappa ,F(\kappa \left)\right)=0
(4)
for x_{
i
} to get the set {\stackrel{\u0304}{x}}_{i}\left(\kappa \right) of fixed point coordinates of the variable x_{
i
} 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 circuitcharacteristic. 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 v_{1} in \stackrel{~}{V}, {f}_{1}({x}_{1}={\kappa}_{1},F({x}_{1}={\kappa}_{1}\left)\right):\mathbb{R}\to \mathbb{R} is a onedimensional characteristic that is called the circuitcharacteristic c(κ_{1}) of the respective SCC. It’s zeros correspond to the fixed point coordinates of variable x_{1}, denoted by \left\{{\stackrel{\u0304}{x}}_{1}\right\}. If we leave the current SCC and go to the next one, we refer to this characteristic as {c}^{k}\left({\kappa}_{1}^{k}\right).

(f)
The corresponding fixed point coordinates of the other variables can be calculated iteratively by inserting the values of the set \left\{{\stackrel{\u0304}{x}}_{1}\right\} into the partial circuitcharacteristics 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 F^{k}.
The structure of the CBA is illustrated in Figure 2 with a flow chart.
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]:
\begin{array}{lll}{\stackrel{\u0307}{x}}_{1}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}& {e}_{N}{x}_{1}& \phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}{r}_{1}{x}_{1}\\ {\stackrel{\u0307}{x}}_{2}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}& \frac{5{x}_{1}}{1+{x}_{1}}\frac{1}{1+{x}_{3}^{4}}{x}_{2}& \phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}{r}_{2}({x}_{1},{x}_{3}){x}_{2}\\ {\stackrel{\u0307}{x}}_{3}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}& \frac{5{x}_{4}}{1+{x}_{4}}\frac{1}{1+{x}_{2}^{4}}{x}_{3}& \phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}{r}_{3}({x}_{2},{x}_{4}){x}_{3}\\ {\stackrel{\u0307}{x}}_{4}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}& \frac{{e}_{M}}{1+{x}_{2}^{4}}{x}_{4}& \phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}{r}_{4}\left({x}_{2}\right){x}_{4}\\ {\stackrel{\u0307}{x}}_{5}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}& \left(\frac{{x}_{1}{x}_{4}}{1+{x}_{1}{x}_{4}}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}\frac{4{x}_{3}}{1+{x}_{3}}\right)\frac{1}{1+{x}_{2}^{4}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{x}_{5}& \phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}{r}_{5}({x}_{1},{x}_{2},{x}_{3},{x}_{4})\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{x}_{5}\\ {\stackrel{\u0307}{x}}_{6}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}& \left(\frac{{x}_{1}{x}_{4}}{1+{x}_{1}{x}_{4}}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}\frac{4{x}_{2}}{1+{x}_{2}}\right)\frac{1}{1+{x}_{3}^{4}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{x}_{6}& \phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}{r}_{6}({x}_{1},{x}_{2},{x}_{3},{x}_{4})\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{x}_{6}\end{array}
(5)
This model describes the interplay between the two lineagespecific counteracting suppressors Gfi1 (x_{2}) and Egr(1,2) (x_{3}) during cellular differentiation for the neutrophil and macrophage cell fate choices, respectively. These are activated by their transcription factors C/EBPα(x_{1}) and PU.1 (x_{4}), respectively. Together, they regulate the expression of lineagespecific downstream genes, which are not further specified in the model and denoted by Mac (x_{5}) and Neut (x_{6}). 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 nondimensionalized after some simplifications by rescaling time and protein concentrations. The two parameters that are left, e_{
N
} and e_{
M
}, 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 μ=e_{
M
} and condition e_{
N
}=e_{
M
} that was created using xpppaut. For e_{
M
}=0 the system has a globally stable fixed point at x=0. The system undergoes a saddlenode bifurcation at {e}_{M}^{\ast}\approx 0.3221. It has a globally stable fixed point for {e}_{M}<{e}_{M}^{\ast} and two stable fixed points divided by an intermediate unstable one for {e}_{M}>{e}_{M}^{\ast}. It can also be seen that the stable fixed point branch that exists for all e_{
M
}represents the neutrophil state, since the fixed point coordinates of the neutrophil specific proteins (x_{1},x_{2},x_{6}) increase monotonically along this branch. The macrophage state is represented by the stable fixed point branch that appears at {e}_{M}^{\ast}.
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 Igraph of system (5) consists of four strongly connected components given by {V}^{1}=\left\{{x}_{1}\right\},{V}^{2}=\{{x}_{2},{x}_{3},{x}_{4}\},{V}^{3}=\left\{{x}_{5}\right\} and {V}^{4}=\left\{{x}_{6}\right\} with circuit sets {C}^{1}={C}^{3}={C}^{4}=\varnothing ,{C}^{2}=\left\{\{{x}_{2},{x}_{4},{x}_{3}\},\{{x}_{2},{x}_{3}\}\right\}, and minimal circuitcovering vertex sets {\stackrel{~}{V}}^{1}={\stackrel{~}{V}}^{3}={\stackrel{~}{V}}^{4}=\varnothing and {\stackrel{~}{V}}^{2}=\left\{{x}_{2}\right\}.
We start with G^{1}(V^{1},E^{1}), which does not contain any circuits. Thus, we just have to solve {\stackrel{\u0307}{x}}_{1}=0 in system (5), which leads to the set {F}^{1}={\stackrel{\u0304}{x}}_{1}^{1}=\left\{{r}_{1}\right\} of fixed points of G^{1}. The fixed point set {}_{u={r}_{1}}{F}^{2}=\left\{{\stackrel{\u0304}{x}}_{2}^{2}\right({x}_{1}={r}_{1}),{\stackrel{\u0304}{x}}_{3}^{2}({x}_{1}={r}_{1}),{\stackrel{\u0304}{x}}_{4}^{2}({x}_{1}={r}_{1}\left)\right\} of G^{2}(V^{2},E^{2}) is calculated by breaking the two circuits at x_{2}, i.e. setting {x}_{2}=:{\kappa}_{1}^{2}. Inserting {\stackrel{\u0304}{x}}_{4}\left({\kappa}_{1}^{2}\right)={r}_{4}\left({\kappa}_{1}^{2}\right) and {\stackrel{\u0304}{x}}_{3}\left({\kappa}_{1}^{2}\right)={r}_{3}({\kappa}_{1}^{2},{r}_{4}({\kappa}_{1}^{2}\left)\right) into {\stackrel{\u0307}{x}}_{2} leads to the circuitcharacteristic
u={r}_{1}{c}^{2}\left({\kappa}_{1}^{2}\right)={r}_{2}\left({r}_{1},{r}_{3}\left({\kappa}_{1}^{2},{r}_{4}\left({\kappa}_{1}^{2}\right)\right)\right){\kappa}_{1}^{2},
(6)
which can, by inserting the respective terms for the synthesis rates r, be rewritten as
{c}^{2}\left({\kappa}_{1}^{2}\right)=\frac{5{r}_{1}}{1+{r}_{1}}\frac{1}{1+{\stackrel{\u0304}{x}}_{3}\left({\kappa}_{1}^{2}\right)}{\kappa}_{1}^{2}
(7)
with
\begin{array}{cc}{\stackrel{\u0304}{x}}_{3}\left({\kappa}_{1}^{2}\right)& =\frac{5{\stackrel{\u0304}{x}}_{4}\left({\kappa}_{1}^{2}\right)}{1+{\stackrel{\u0304}{x}}_{4}\left({\kappa}_{1}^{2}\right)}\frac{1}{1+{\left({\kappa}_{1}^{2}\right)}^{4}}\phantom{\rule{1em}{0ex}}\text{and}\\ {\stackrel{\u0304}{x}}_{4}\left({\kappa}_{1}^{2}\right)& =\frac{{e}_{M}}{1+{\left({\kappa}_{1}^{2}\right)}^{4}}.\end{array}
(8)
This characteristic is shown in Figure 5 (center row), along with the sets {\stackrel{\u0304}{x}}_{i}\left({\kappa}_{1}^{2}\right), i=3,4,5,6, for parameter values e_{
M
}={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 x_{3},x_{4},x_{5}and x_{6}behave monotonically with the input κ_{2}, which represents the neutrophil state. The macrophage specific proteins x_{3},x_{4}and x_{5}decrease with increasing {\kappa}_{1}^{2}, x_{6}increases.

2.
Looking at the characteristics (center row) for different values e_{
M
}, it is monotonically decreasing for {e}_{M}<{e}_{M}^{\ast} (left), and thus has a single zero, which corresponds to the single fixed point branch for {e}_{M}<{e}_{M}^{\ast}. For the value e_{
M
}=0.2, which is chosen here, we get the fixed point \stackrel{\u0304}{x}(\mu =0.2)={0.2,0.81,0.43,0.14,0.86,1.75}, as indicated in the graphs. This state represents an intermediate nondifferentiated progenitor cell state. The saddlenode bifurcation is represented by the second zero of the characteristic that appears at {e}_{M}={e}_{M}^{\ast} (center column). The respective fixed point set is {\stackrel{\u0304}{x}}^{1}(\mu =0.3221)=\{0.3221,0.51,1.09,0.30,2.04,0.59} and {\stackrel{\u0304}{x}}^{2}(\mu =0.3221)=\{0.3221,1.21,0.15,0.10,0.17,2.22}.Finally, the characteristic has three zeros for {e}_{M}>{e}_{M}^{\ast} (right column) and thus the system has three fixed points in this range. For the chosen value e_{
M
}=0.5 we can read off the fixed point set {\stackrel{\u0304}{x}}^{1}(\mu =0.5)=\{0.5,0.20,1.67,0.5,2.7,0.03\}, {\stackrel{\u0304}{x}}^{2}(\mu =0.5)=\{0.5,0.75,1.05,0.38,1.68,0.84\} and {\stackrel{\u0304}{x}}^{3}(\mu =0.5)=\{0.5,1.66,0.03,0.06,0.02,2.52\}. Here, {\stackrel{\u0304}{x}}^{1} represents the macrophage state, where Egr and PU.1 are highly expressed, and C/EBPαis low, {\stackrel{\u0304}{x}}^{3} stands for the neutrophil state in which C/EBPαis dominant, and {\stackrel{\u0304}{x}}^{2} is an unstable intermediate state that separates the two basins of attraction.
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 xaxis without intersection are bifurcation value candidates. Furthermore, looking at this example, a selfevident guess would be to assume that stability can be determined in the same way as for onedimensional vector fields: The fixed points are stable if the slope of the characteristic at the corresponding zero κ^{∗} is negative, i.e. \mathrm{dc}\left(\kappa \right)/\mathrm{d\kappa}{\left\right.}_{\kappa ={\kappa}^{\ast}}<0, and it is unstable if the slope is positive, i.e. \mathrm{dc}\left(\kappa \right)/\mathrm{d\kappa}{\left\right.}_{\kappa ={\kappa}^{\ast}}>0. We will further investigate these assumptions in the following subsections. In order to do so, we consider in the following strongly connected Igraphs, 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 {\kappa}_{1}\in \mathbb{R} is the value of the variable x_{1}, 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 J_{
f
}(x) of the system:
Proposition 1
\frac{\mathrm{dc}\left({\kappa}_{1}\right)}{d{\kappa}_{1}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\frac{d{f}_{1}\left(x\right)}{d{x}_{1}}{\left\right.}_{({x}_{1},F({x}_{1}\left)\right)}
(9)
\phantom{\rule{0.5em}{0ex}}=\phantom{\rule{0.5em}{0ex}}det\phantom{\rule{0.5em}{0ex}}{J}_{f}\left(x\right){\left\right.}_{({x}_{1},F({x}_{1}\left)\right)}\xb7{det}^{1}{J}_{f}^{V\setminus \left\{{v}_{1}\right\}}\left(x\right){\left\right.}_{F\left({x}_{1}\right)}
(10)
with F(x_{1}) denoting the fixed point coordinates of variables x_{2},…,x_{V}in dependence of x_{1}, and {J}_{f}^{V\setminus \left\{{v}_{1}\right\}}\left(x\right) is the Jacobian matrix of the subnetwork spanned by the vertices V∖{v_{1}}.
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 {\kappa}_{1}^{\ast} with c\left({\kappa}_{1}^{\ast}\right)=0, and we will in the following subsection sometimes denote the corresponding fixed point with \stackrel{\u0304}{x}, if appropriate.
Instability of fixed points
From Proposition 1 it follows that a positive slope \frac{\mathrm{dc}\left({\kappa}_{1}\right)}{d{\kappa}_{1}}{\left\right.}_{{\kappa}_{1}^{\ast}}>0 implies that det\phantom{\rule{0.3em}{0ex}}{J}_{f}\left(\stackrel{\u0304}{x}\right){\left\right.}_{({\stackrel{\u0304}{x}}_{1},F({\stackrel{\u0304}{x}}_{1}\left)\right))} and det\phantom{\rule{0.3em}{0ex}}\underset{f}{\overset{V\setminus \left\{{v}_{1}\right\}}{J}}\left(\stackrel{\u0304}{x}\right){\left\right.}_{F\left({\stackrel{\u0304}{x}}_{1}\right)} have the same signs. According to the HartmanGrobman theorem (see e. g. [17]), a fixed point \stackrel{\u0304}{x} is unstable if {J}_{f}\left(\stackrel{\u0304}{x}\right) has at least one eigenvalue with positive real part. Unfortunately, we are not aware of a relation between the determinant of J_{
f
}(x) and it’s minors that can be used to show the following: If det\phantom{\rule{0.3em}{0ex}}{J}_{f}\left(\stackrel{\u0304}{x}\right) and det\phantom{\rule{0.3em}{0ex}}\underset{f}{\overset{V\setminus \left\{{v}_{1}\right\}}{J}}\left(\stackrel{\u0304}{x}\right) have the same signs, this implies the existence of an eigenvalue with positive real part and hence implies instability of \stackrel{\u0304}{x}. 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 \stackrel{~}{V} that consists of one single element v_{1}. 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 circuitclosing step. The Igraph 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 circuitclosing 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 circuitcharacteristic 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:
\begin{array}{ll}{\stackrel{\u0307}{m}}_{i}& =\frac{\alpha}{1+{p}_{j}^{n}}+{\alpha}_{0}{m}_{i}=:r\left({p}_{j}\right){m}_{i}\phantom{\rule{2em}{0ex}}\\ {\stackrel{\u0307}{p}}_{i}& =\beta ({m}_{i}{p}_{i}),\phantom{\rule{2em}{0ex}}\end{array}
(11)
with i={lacI,tetR,cI}, j={cI,lacI,tetR}, and m_{
i
} and p_{
i
} 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 {m}_{i}^{\mathrm{min}}={\alpha}_{0}{m}_{i}^{\mathrm{max}}=\alpha +{\alpha}_{0} and {p}_{i}^{\mathrm{min}/\mathrm{max}}={m}_{i}^{\mathrm{min}/\mathrm{max}}i=1,2,3. The Igraph (Figure 6) is strongly connected, the circuit set C consists of one subset that contains all six nodes, C={{m_{
i
}p_{
i
}}_{i=1,2,3}}, and hence the set \stackrel{~}{V} has one single element and the graph is a LVG.
Note that because of the symmetry of the model, the circuitcharacteristic is independent of the choice of \stackrel{~}{V} here. It is given by
c\left({\kappa}_{1}\right)=r\left({\stackrel{\u0304}{p}}_{3}\right({\stackrel{\u0304}{m}}_{3}\left({\stackrel{\u0304}{p}}_{2}\right({\stackrel{\u0304}{m}}_{2}\left({\stackrel{\u0304}{p}}_{1}\right({\kappa}_{1}\left)\right)\left)\right)\left)\right){\kappa}_{1},
(12)
which can be simplified to
c\left({\kappa}_{1}\right)=r\left(r\right(r\left({\kappa}_{1}\right)\left)\right){\kappa}_{1},
(13)
where we have used {\stackrel{\u0304}{p}}_{i}\left({\stackrel{\u0304}{m}}_{i}\right)={\stackrel{\u0304}{m}}_{i}, {\stackrel{\u0304}{m}}_{i}\left({\stackrel{\u0304}{p}}_{j}\right)=r\left({\stackrel{\u0304}{p}}_{j}\right) and {\stackrel{\u0304}{p}}_{1}={\stackrel{\u0304}{p}}_{2}={\stackrel{\u0304}{p}}_{3}. 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 {J}_{f}\left(\stackrel{\u0304}{x}\right) are plotted as functions of β for parameter values α=290, n=2 and α_{0}=10. For these parameter values the system has a fixed point {\stackrel{\u0304}{m}}_{i}={\stackrel{\u0304}{p}}_{i}=12 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 {J}_{f}\left(\stackrel{\u0304}{x}\right) and in particular the stability of \stackrel{\u0304}{x} 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.
As a further example we consider the tryptophan regulation network in Escherichia coli described in [11], which can be written as
\phantom{\rule{16.0pt}{0ex}}\begin{array}{lll}{\stackrel{\u0307}{x}}_{1}\phantom{\rule{0.3em}{0ex}}=& {k}_{1}{O}_{t}C({x}_{4},\phantom{\rule{0.3em}{0ex}}{t}_{1},{m}_{1})\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}({\gamma}_{1}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}\mu ){x}_{1}& \phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}{r}_{1}\left(x\right)\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}({\gamma}_{1}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}\mu ){x}_{1}\\ {\stackrel{\u0307}{x}}_{2}\phantom{\rule{0.3em}{0ex}}=& {k}_{2}{x}_{1}C({x}_{4},\phantom{\rule{0.3em}{0ex}}{t}_{2},\phantom{\rule{0.3em}{0ex}}{m}_{2})\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}({\gamma}_{2}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}\mu ){x}_{2}& ={r}_{2}\left(x\right)\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}({\gamma}_{2}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}\mu ){x}_{2}\\ {\stackrel{\u0307}{x}}_{3}=& {k}_{3}{x}_{2}\mu {x}_{3}& ={r}_{3}\left(x\right)\mu {x}_{3}\\ {\stackrel{\u0307}{x}}_{4}=& {k}_{4}C({x}_{4},{t}_{3},{m}_{3}){x}_{3}g\frac{{x}_{4}}{{x}_{4}+{K}_{g}}\mu {x}_{4}& ={r}_{4}\left(x\right)\mu {x}_{4},\end{array}
(14)
where the state vector x corresponds to the free operator sites (O_{
R
}), mRNA (M), enzyme (E) and tryptophan (T) concentrations. C(x K m) are sigmoidally decreasing functions,
C(x,K,m)=\frac{{K}^{m}}{{K}^{m}+{x}^{m}}.
(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(x_{4}t_{1}m_{1}), mRNA attenuation (C(x_{4}t_{2}m_{2})) and enzyme inhibition (C(x_{4}t_{3}_{m3})). The parameters k_{1}k_{2}k_{3}and k_{4}represent kinetic rate constants for synthesis of free operator, mRNA transcription, translation and tryptophan synthesis, respcetively, K are halfsaturation constants for the inhibition processes, O_{
t
} denotes the total operator site concentration, and γand μrefer to degradation and diluation rates due to cell growth. The term g\frac{{x}_{4}}{{x}_{4}+{K}_{g}} describes the uptake of tryptophan for protein synthesis in the cell.
Analyzing this system with the parameter values given in [11] ({k}_{1}=50{\mathrm{min}}^{1},{O}_{t}=3.32\mathrm{nM},{t}_{1}=3.53\mathrm{\mu M},{m}_{1}=1.92,{\gamma}_{1}=0.5{\mathrm{min}}^{1},\mu =0.01\mathrm{mi}{n}^{1},{k}_{2}=15\mathrm{mi}{n}^{1},{t}_{2}=0.04\mathrm{\mu M},{m}_{2}=1.72,{\gamma}_{2}=15{\mathrm{min}}^{1},{k}_{3}=90\mathrm{mi}{n}^{1},{k}_{4}=59{\mathrm{min}}^{1},{t}_{3}=810\mathrm{\mu M},{m}_{3}=1.2,g=25μM K_{
g
}=0.2μM) using xppaut and choosing the dilution rate μas bifurcation parameter, the system shows a Hopf bubble between {\mu}_{1}^{\ast}=0.02486 and {\mu}_{2}^{\ast}=0.1529 (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.
The corresponding Igraph is shown in Figure 9. It is strongly connected.
The circuit set C and the minimal circuit covering vertex set \stackrel{~}{V} are C=\left\{\left\{{x}_{4}\right\},\{{x}_{2},{x}_{3},{x}_{4}\},\{{x}_{1},{x}_{2},{x}_{3},{x}_{4}\}\right\} and \stackrel{~}{V}=\left\{{x}_{4}\right\}. Since \stackrel{~}{V} consists of a single element, this system is a LVG and only one circuitclosing step is necessary to calculate the set of fixed points. The circuitcharacteristic can be calculated analytically here and is given by
c\left({\kappa}_{4}\right)={r}_{4}\left({\stackrel{\u0304}{x}}_{3}\right({\kappa}_{4}),{\kappa}_{4})\mu {\kappa}_{4},
(16)
where r_{4}can iteratively be calculated via
\begin{array}{ll}{\stackrel{\u0304}{x}}_{1}\left({\kappa}_{4}\right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}& =\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\frac{{r}_{1}\left({\kappa}_{4}\right)}{{\gamma}_{1}+\mu}\phantom{\rule{2em}{0ex}}\\ {\stackrel{\u0304}{x}}_{2}\left({\kappa}_{4}\right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}& =\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\frac{{r}_{2}\left({\stackrel{\u0304}{x}}_{1}\right({\kappa}_{4}),{\kappa}_{4})}{{\gamma}_{2}+\mu}\phantom{\rule{2em}{0ex}}\\ {\stackrel{\u0304}{x}}_{3}\left({\kappa}_{4}\right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}& =\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\frac{{r}_{3}\left({\stackrel{\u0304}{x}}_{2}\right({\kappa}_{4}\left)\right)}{\mu}.\phantom{\rule{2em}{0ex}}\end{array}
(17)
As can be seen in Figure 10, c(κ_{4}) is strictly decreasing (bottom row), since all circuits in the graph are negative.
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 {\kappa}_{4}^{\ast}=31.8, which corresponds to the fixed point coordinate {\stackrel{\u0304}{x}}_{4}. Inserting this value into {\stackrel{\u0304}{x}}_{i}\left({\kappa}_{4}\right), i=1,2,3, we get the fixed point \stackrel{\u0304}{x}(\mu =0.01)=(4.71,4.82\xb71{0}^{5},0.43,31.82), and likewise for the other dilution rates. The qualitative courses of c(κ_{4}) and also for the fixed point sets {\stackrel{\u0304}{x}}_{i}\left({\kappa}_{4}\right) 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.