 Methodology article
 Open Access
Structurally robust biological networks
 Franco Blanchini^{1} and
 Elisa Franco^{2}Email author
https://doi.org/10.1186/17520509574
© Blanchini and Franco; licensee BioMed Central Ltd. 2011
Received: 10 November 2010
Accepted: 17 May 2011
Published: 17 May 2011
Abstract
Background
The molecular circuitry of living organisms performs remarkably robust regulatory tasks, despite the often intrinsic variability of its components. A large body of research has in fact highlighted that robustness is often a structural property of biological systems. However, there are few systematic methods to mathematically model and describe structural robustness. With a few exceptions, numerical studies are often the preferred approach to this type of investigation.
Results
In this paper, we propose a framework to analyze robust stability of equilibria in biological networks. We employ Lyapunov and invariant sets theory, focusing on the structure of ordinary differential equation models. Without resorting to extensive numerical simulations, often necessary to explore the behavior of a model in its parameter space, we provide rigorous proofs of robust stability of known biomolecular networks. Our results are in line with existing literature.
Conclusions
The impact of our results is twofold: on the one hand, we highlight that classical and simple control theory methods are extremely useful to characterize the behavior of biological networks analytically. On the other hand, we are able to demonstrate that some biological networks are robust thanks to their structure and some qualitative properties of the interactions, regardless of the specific values of their parameters.
Keywords
 Lyapunov Function
 Biological Network
 Robust Stability
 MAPK Cascade
 Inducer Concentration
Background
The complex biochemistry of living organisms very often outperforms electrical and mechanical devices in terms of adaptability and robustness. Mapping such intricate reaction networks to high level design principles is the goal of systems biology, and it requires an immense collaborative effort among different disciplines, such as physics, mathematics and engineering [1].
The most classical example of robust molecular circuitry is probably given by bacterial chemotaxis [2, 3]. The action of the flagellar motor of E. coli cells is driven by a cascade of signaling proteins, whose active or inactive state is determined by the presence of nutrient in the environment. Both analysis on a simplified ordinary differential equation (ODE) model [2] and experiments [3] showed how the flagellar motion of E. coli presents a robustly stable steady state: steps in the nutrient concentration only temporarily alter the motor equilibrium. Cells are therefore sensitive to nutrient gradients, but always return to their stable motion mode (such property is also referred to as adaptability). Such stable steady state is only a function of the concentrations of the signaling cascade protein components and a few binding rates, and is therefore independent of external inputs. Further analysis also demonstrated how integral feedback is present in the chemotaxis network, and guarantees robustness (perfect adaptation) of the equilibrium [4].
In this work, we are going to ask a simple question: are there biological systems that present structurally stable equilibria, and preserve this property robustly with respect to their specific parameters? This question has been considered before in the literature. For instance, through extensive numerical analysis on threenode networks, the authors of [5] have shown that adaptability of these systems can be investigated solely based on their structure, regardless of the chosen reaction parameters. In [6], through numerical exploration of the Jacobian eigenvalues for two, three and four node networks, the authors isolated a series of interconnections which are stable, robustly with respect to the specific parameters. Such structures also turned out to be the most frequent topologies in existing biological networks databases. Numerical simulation has arguably been the most popular tool to investigate robustness of biological networks [7–12]. Analytical approaches to the study of robustness have been proposed in specific contexts. A series of recent papers [13, 14] focused on input/output robustness of ODE models for phosphorylation cascades. In particular, the theory of chemical reaction networks is used in [14] as a powerful tool to demonstrate the property of absolute concentration robustness. Indeed, the socalled deficiency theorems [15] are to date some of the most general results to establish robust stability of a chemical reaction network. Monotonicity is also a structural property that is useful to demonstrate robust dynamic behaviors of a class of biological models [16, 17]. Robustness has also been investigated in the context of compartmental models, which are often encountered in biology and chemistry [18].
In this work, we provide a simple and general theoretical tool kit for the analysis biomolecular systems. Such tools are suitable for the investigation of robust stability by means of Lyapunov and setinvariance methods. Provided that certain standard properties are verified, we demonstrate how a number of well known biological networks are asymptotically stable, robustly with respect to the model parameters. In some cases, we are also able to provide robust bounds on the system performance. Our approach does not require numerical simulation efforts. The contribution of the paper can be summarized as follows.

The framework we suggest is easy and intuitive for biologists to formulate qualitative models without the need of exact mathematical expressions and parameters. We will propose analytical methods that only rely on qualitative interactions between network components.

The properties that can be derived from such modeling are, consequently, structurally robust because they are not inferred from mathematical formulas arbitrarily chosen to fit data.

We suggest techniques based on setinvariance and Lyapunov theory, in particular piecewiselinear functions, to show that such models are amenable for robust investigation by engineers and mathematicians. Such techniques are believed to be quite effective and promising in dealing with biological robustness [19], [20].

We consider several models from the literature, reporting the original equations, and rephrasing them in our setup as case studies.

We show how robust certifications can be given to important properties (some of which have been established based on specific models).
Methods
Robustness
where x is the system state, u models external inputs, and both are vectors of appropriate dimensions. Such class of models is appropriate for biological systems where stochasticity and anisotropy can be neglected. We define robustness as follows:
Definition 1 Let be a class of systems and be a property pertaining such a class. Given a family we say that is robustly verified by ℱ, in short robust, if it is satisfied by each element of ℱ.
with positive and constant coefficients a, b, c, d. Assume . Then we can say that is robust. The situation is different if we admit that a(t), b(t), c(t), d(t) can vary with time, yielding a system which is possibly unstable.
where is fixed and κ (·) is a negative definite function (i.e. κ (·) < 0 on all its domain, except for κ (0) = 0).
Nonsmooth Lyapunov functions
Positively invariant sets
We are interested in cases where the trajectories of system (1) remain trapped in bounded sets at all times, therefore behaving consistently with respect to some desired criterion.
We say that a subset of the state space is positively invariant if implies that also for all t > 0. The following theorem (which relies on the concept of Lyapunov function) provides a general necessary and sufficient condition for a set to be invariant.
For instance, if our constraining functions are linear, s^{⊤}x ≤ σ, the Nagumo conditions are . We refer the reader to [22] for further details on positively invariant sets; more recent works on this topic are [23] and [24].
Structural robustness investigation for biological networks
RNA species x_{2} determines the production rate of protein x_{3} by indexing the corresponding reaction rate as a_{32}. Following the standard notation in control theory, we assume that the production rate of protein x_{1} is driven by some external signal or input u_{1}, and that RNA u_{2} also acts as an external input on RNA x_{2}. We assume that all the system parameters are positive and bounded scalars. Terms a_{ ij } are first order production rates: species i is produced at a rate which is linear in species j; b_{ ih } denote in this case first order degradation rates. The term d_{21}(x_{1}) is a well known Hill function term [26]. The stability properties of this small network can be immediately assessed: x_{1} will converge to its equilibrium . Similarly, , . Regardless of the specific parameter values, and therefore robustly, the system is stable. The equilibrium could grow unbounded with u_{1}, however is always bounded.
We remark that the knowledge of functions a_{ ij }x, b_{ ih }x and d(·) is not necessary at all: the previous conclusions can be easily derived by the qualitative information that d_{21} is strictly decreasing and asymptotically converging to 0, while b_{11}x_{1}, b_{22}x_{2}, , a_{32}x_{2} and b_{33}x_{3} are increasing.
It is appropriate at this point to outline a series of general assumptions that will be useful in the following analysis.
where x_{ i } , i = 1,..., n are the dynamic variables. For the sake of notation simplicity, we are not denoting external inputs with a different symbol. Inputs can be easily included as dynamic variables which are not affected by other states and have the desired dynamics. The sets denote the subsets of variables affecting x_{ i } . The different terms in equation (3) are associated with a specific biological and physical meaning. The terms f_{ ij } (· , ·) are associated with production rates of reagents, typically, these functions are assumed polynomial in their arguments; similarly, terms g_{ ih } (· , ·) model degradation or conversion rates and are also likely to be polynomial in practical cases. Finally, terms c(·) and d(·) are associated with monotonic nonlinear terms, often given by MichaelisMenten or Hill functions [26]. We assume that system (3) satisfies the following assumptions:
A1 (Smoothness) Functions f_{ ij } (· , ·), g_{ ih } (· , ·), c_{ is } (·) and d_{ il } (·) are unknown nonnegative continuously differentiable functions.
A2 f_{ ij } (x_{ i } , 0) = 0 and g_{ ih } (x_{ i } , 0) = 0, ∀x.
A4 (Saturation) Functions c_{ is } (x_{ s } ) and d_{ il } (x_{ l } ) are nonnegative and, respectively, nondecreasing and nonincreasing. Moreover c_{ is } (∞) > 0 and d_{ il } (0) > 0.
A5 Functions g_{ ih } (· , ·) are null at the lower saturation levels : g_{ ih } (0, x_{ h } ) = 0, ∀x_{ h } .
In view of the nonnegativity assumptions and Assumption 5, the general model (4) is a nonlinear positive system, according to the next proposition, and its investigation will be restricted to the positive orthant.
Proposition 1 The nonnegative orthant x_{ i }≥ 0 is positively invariant for system (4).
The above expression is always valid due to the smoothness assumption A1 (see [18], Section 2.1).
To simplify the notation, we have considered functions depending on two variables at most. However, we can straightforwardly extend assumptions A1A5 to multivariate functional terms in equation (3). In turn, the model structure (4) can be easily generalized to include terms as a(x_{ i } , x_{ j } , x_{ k } ,...), b(x_{ i } , x_{ j } , x_{ k } ,...), c(x_{ i } , x_{ j } , x_{ k } ,...), d(x_{ i } , x_{ j } , x_{ k } ,...).
If we restrict our attention to the general class of models (4), under assumptions A1A5, we can proceed to successfully analyzing the robust stability properties of several biological network examples.
The structural analysis of system (4) can be greatly facilitated whenever it is legitimate to assume that functions a, b, c d have certain properties. For the reader's convenience, a list of possible properties is given below. Given a general function f(x):
P1 f (x) = const ≥ 0 is nonnegativeconstant.
P2 f(x) = const > 0 is positiveconstant.
P3 f (x) is sigmoidal: it is nondecreasing, f(0) = f '(0) = 0, if 0 < f(∞) < ∞ and its derivative has a unique maximum point, for some .
P4 f (x) is complementary sigmoidal: it is nonincreasing, 0 < f(0), f'(0) = 0, f(∞) = 0 and its derivative has a unique minimum point. In simple words, f is a CSM function iff f(0)  f(x) is a sigmoidal function.
P5 f (x) is constantsigmoidal, the sum of a sigmoid and a positive constant.
P6 f (x) is constantcomplementarysigmoidal, the sum of a complementary sigmoid and a constant.
P7 f (x) is increasingasymptoticallyconstant: f'(x) > 0, 0 < f (∞) < ∞ and its derivative is decreasing.
P8 f (x) is decreasingasymptoticallynull: f'(x) < 0, f (∞) = 0 and its derivative is increasing.
P9 f (x) is decreasingexactlynull: f'(x) < 0, for and f(x) = 0 for for some .
P10 f (x) is increasingasymptoticallyunbounded: f'(x) > 0, f (∞) = + ∞.
Network graphs
Building a dynamical model for a biological system is often a long and challenging process. For instance, to reveal dynamic interactions among a pool of genes of interest, biologists may need to selectively knockout genes, set up micro RNA assays, or integrate fluorescent reporters in the genome. The data derived from such experiments are often noisy and uncertain, which implies that also the estimated model parameters will be uncertain. However, qualitative trends can be reliably assessed in the dynamic or steady state correlation of biological quantities.
Graphical representations of such qualitative trends are often used by biologists, to provide intuition regarding the network main features. We believe that, indeed, such graphs may be useful even to immediately construct models analogous to (3). We propose a specific method to construct such graphs: the biochemical species of the network will be associated to the nodes in the graph, the qualitative relationships between the species will be instead associated with different types of arcs: in particular, the terms of a, b, c and d will be represented as arcs having different endarrows, as shown in Figure 1A. These graphs can be immediately constructed, by knowing the correlation trends among the species of the network, and aid the construction of a dynamical model. For simple networks, this type of graph may provide intuition regarding their behavior and may facilitate their structural robustness analysis. For instance, the graph associated to equations (2) is shown in Figure 1B. Throughout the paper, we will consider similar case studies and use their graph representation as a visual support for the analysis.
Remark 1 In this work, properties such as positivity, monotonicity, boundedness and other functional characteristics are labeled as "qualitative and structural properties"[27]. Through such properties, we can draw conclusions on the dynamic behaviors of the considered systems without requiring specific knowledge of parameters and without numerical simulations. However, it is clear that our approach requires more information than other methods, such as boolean networks and other graphbased frameworks.
Investigation method
The main objective of this work is to show that, at least for reasonably simple networks, structural robust stability can be investigated with simple analytical methods, without the need for extensive numerical analysis. We will suggest a two stage approach:

Preliminary screening: establish essential information on the network structure, recognizing which properties (such as P1P10) pertain to each link.

Analytical investigation: infer robustness properties based on dynamical systems tools such as Lyapunov theory, set invariance and linearization.
Results and Discussion
In this section we will analyze five biological networks as case studies. Three of such examples, the Larabinose, the sRNA and the Lac Operon networks, model the interaction and control of expression of a set of genes. The cAMP and the MAPK pathways are instead signaling networks, namely they represent sets of chemical species interacting for transmission and processing of upstream input signals. These networks are all wellknown in the literature, and have been characterized mainly through experimental and numerical methods, although the MAPK pathway, for instance, has been thoroughly analyzed using the theory of monotone systems [17].
We will provide rigorous proofs that these networks are either mono or multistable in a robust manner. Such demonstrations rely on Lyapunov functions and invariant sets theory, according to our proposed methodology. In some cases, we are also able to provide bounds on their speed of convergence.
The Larabinose network
where u is nonnegativeconstant, c_{1}, b_{11} and b_{22} are positiveconstant, while c_{1u}(u) and c_{2u 1}(u) are sigmoidal with respect to u, the latter increasing with respect to x_{1}. The graph representation of this network is in Figure 3A.
For this elementary network the analysis is straightforward. Variable x_{1} is not affected by x_{2}. Since c_{1u}(u) is bounded, x_{1} is also bounded and converges to an equilibrium point which is monotonically increasing in u. In turn, x_{2} is also positive and bounded for any value of u and stably converges to a unique equilibrium point , which is a monotonically increasing function of u (partially activated by ). The positive term c_{1} prevents x_{1}(t) and x_{2}(t) from staying at zero. It is worth remarking that the hierarchical structure of this network greatly facilitates the analysis; equilibria can in fact be iteratively found and their stability properties characterized.
The sRNA pathway
thus the function represents the worst case between the sum and the mismatch.
The following proposition shows that the sRNA pathway is a typical system in which robustness is structurally assured. We report the full demonstration of this proposition, because its steps and the techniques used are a model for the subsequent proofs in this paper.
for some β > 0 and any x_{1}(0) ≥ 0, x_{2}(0) ≥ 0. Moreover, no oscillations are possible around the equilibrium, in the sense that the condition or occurs at most once.
Define κ = (c_{1} + c_{2})/min {b_{11}, b_{22}} then for s(x_{1}, x_{2}) > κ the derivative becomes negative so s(x_{1}, x_{2}) cannot exceed κ (See Theorem 1).
where we have subtracted the null terms (10) and where we have exploited the fact that b_{12}(x_{1}, x_{2})x_{1} = b_{21}(x_{1}, x_{2})x_{2} is increasing in both variables. The inequality (CPD in Figure 2B) can be similarly proved to hold in the sector and .
with β = min{b_{11}, b_{22}}. This implies (9) and the uniqueness of the equilibrium point.
We finally need to show that there are no oscillations. To this aim we notice that the sectors DPA, and , and its opposite
CPB, and , are both positive invariant sets.
hence the claimed invariance of sector DPA. The proof of the invariance of sector CPB is identical.
Remark 2 Note that the constructed Lyapunov function does not depend on the system parameters. This fact can be used to prove that if the transcription rates c_{1}(t) and c_{2}(t) are timevarying, but bounded, we have convergence to a neighborhood whose amplitude, obviously, depends on the bounds of c_{1}(t) and c_{2}(t). It is realistic to assume that the transcription rates vary over time: for instance, if the environmental conditions change, the cell may need to down or upregulate entire groups of transcripts and therefore increase or decrease c_{2}(t).
The following corollary evidences the positive influence of c_{2}, which is positive over x_{2} and negative over x_{1}.
Corollary 1 Assume that x_{1}(0), x_{2}(0) is at the steady state corresponding to and . Consider the new input (keeping ). Then the system converges to a new equilibrium with and . There is no undershoot, respectively, overshoot.
where (keep in mind that b_{21}(x_{1}, x_{2})x_{1} = b_{12}(x_{1}, x_{2})x_{2}). The sign of the entries in the second column are negative and positive respectively, therefore, the steadystate values and are decreasing and increasing functions of c_{2}.
The absence of overshoot and undershoot is an immediate consequence of the invariance of the sector and previously proved.
Obviously, decreasing c_{2} increases x_{1} and decreases x_{2} and the same holds if we commute 1 and 2. It is worth noting that the same conclusions regarding the lack of multistability and oscillations for the sRNA pathway may be reached by qualitative analysis of the system's nullclines.
The cAMP dependent pathway
Our preliminary analysis allows us to assume: a_{1u}, a_{23}: decreasingexactlynull with threshold values and ; d_{32}, a_{31}: decreasingasymptoticallynull, b_{32} and g_{33} = b_{33}(x_{3})x_{3}: increasingasymptoticallyconstant; b_{11}, b_{22}are positiveconstant.
The solution of the previous equation is uniquely defined for each u since the function ξ^{1}(x_{1}) on the right is strictly increasing and grows to infinity, precisely . Biologically, this means that external glucose signals are mapped to internal active Gprotein concentration with a bijection, before saturating.
Also we have to note that the model is consistent with mass conservation: since a_{1u}(x_{1}) and a_{23}(x_{2}) are zero above the thresholds and , we have that and for and , respectively; therefore we assume , , for all t ≥ 0.
where as previously defined . All the equilibrium values , and are increasing functions of u. If condition (13) is satisfied , the equilibrium is unique and locally stable.
The previous proposition assures only local stability, but this result can be extended to global stability. To this aim, we will assume that x_{1} is at its equilibrium value . Furthermore, under a suitable condition a performance bound on the transient values of x_{3}(t) can be given.
The proof can be found in Section S{II of the Additional File 1.
Remark 3 The condition (14) has the following interpretation. It basically states that the inhibiting term b_{33}(x_{3})x_{3}at "full force" (x_{3}suitably large) dominates the activating term d_{32}(x_{2}) + a_{31}(x_{2})ξ when x_{2}is small. Note that, indeed, the feedback terms modulated by the two phosphodiesterases act in a complementary manner, in order to maintain a bounded concentration of cAMP in the cell.
Remark 4 The system, even if initialized with small values x_{ 2 } (0) and x_{3}(0), may exhibit a spike of cAMP x_{3}which is bounded by (15), if condition (14) is satisfied. If x_{ 3 }(0) is small, then the bound is d_{32}(0) + a_{31}(0)ξ (u): the amplitude of the spike is, in general, an increasing function of the glucose concentration u. If condition (14) fails, then (see Figure S2 in the Additional File) the spike of x_{3}(t) can be arbitrarily large; thus condition (14) can be seen as a threshold.
The Lac operon
This genetic network was originally studied by Monod and Jacob [33]. The natural nutrient for E. coli bacterial cells is glucose, which is metabolized by enzymes normally produced by the bacteria. When glucose is absent, but the allolactose inducer is present in their environment, E. coli activates a set of genes that will regulate the lactose intake and breakdown. In particular, the cells start producing a permease protein, which binds to the cell membrane and increases the inflow of lactose; and cells also start producing the βgalactosidase protein, which converts lactose in allolactose.
In this paper we will consider the deterministic model proposed in [34]. This simple model does not capture the stochasticity of this genetic circuit, but it does explain the bimodal behavior of the system. Such behavior is observable experimentally: within the same population, the operon can be either induced or uninduced. Our analysis shows that for low or high intracellular inducer concentrations, the system is monostable and respectively reaches an uninduced or induced equilibrium; however, at intermediate inducer concentrations the system becomes multistable.
where c_{13}(x_{3}) = f_{1}(x_{3}), b_{11} = δ_{1}, a_{21} = β_{1}, b_{22} = δ_{2}, a_{32}(u) = f_{2}(u) =, b_{32}(x_{3}) = f_{3}(x_{3}), c_{3u}= β_{2}, b_{33} = δ_{3}, c_{43}(x_{3}) = γ f_{1}(x_{3}) and b_{44} = δ_{4}. This corresponds to the network in Figure 3C.
From our preliminary analysis step: c_{13} is constantsigmoidal, a_{32}(u) and b_{32}(x_{3}) are increasingasymptoticallyconstant, and the remaining functions a_{21}, b_{11}, b_{22} and b_{33} are positiveconstant.
We can start to study this network without any specific knowledge of the parameters in equations (17). First of all, as evident in Figure 3C, note that the βgalactosidase concentration x_{4} does not affect any other chemical species: therefore, the fourth equation can be considered separately. As long as the inducer concentration of x_{3} within the cell reaches an equilibrium , x_{4} converges to . Therefore, we can restrict our attention to the first three equations; this is consistent with the model proposed in [35, 36]. From now on we will consider this reduced model (see Section SIII of the Additional File), neglecting the linear term c_{3u}u as in [35, 36]. We will not introduce delays in our model, as done in [37]. Our preliminary screening also shows that the evolution of this system is necessarily bounded. Indeed x_{1} receives a bounded signal from x_{3} and the degradation term b_{11}x_{1} keeps x_{1} bounded. In turn, x_{2} remains bounded. The inducer concentration x_{3} receives a bounded signal form u and x_{2}; therefore x_{3} stays bounded as well, being both a_{32}(u) and b_{32}(x_{3}) bounded.
The following proposition evidences that fundamental results can be established starting from our general framework. These results are consistent with the findings in [36], whose analysis relies on assuming Hilltype functions in the model.
Proposition 5 For any functional terms in Equations 17, satisfying the general assumptions formulated above, the system admits a unique equilibrium for large u > 0 or small u > 0.
For some chioces of such functional terms, the system may have multiple positive equilibria x^{ A } , x^{ B } , x^{ C } ,... ∈ IR^{3}(typically three) for intermediate values of u. If multiple equilibria exist, then they are ordered in the sense that x^{ A } ≤ x^{ B } ≤ x^{ C } ... where the inequality has to be considered componentwise. If the equilibria are all distinct, then they are alternatively stable and unstable. In the case of three equilibria, x^{ A } , x^{ B } , x^{ C } they are stable, unstable and stable, respectively. Finally, given any equilibrium point, the positive and negative cones x ≤ x* and × ≥ x* are positively invariant.
The proof is given in Section SIII of the Additional File. The coneinvariance property implies that the state variables cannot exhibit oscillations around their equilibria. For instance, if x^{ A } is the first (hence stable) equilibrium, given any initial condition upper bounded by x^{ A } (x(0) x^{ A } ) in the domain of attraction, the convergence to x^{ A } has no overshoot (and if x(0) ≥ x^{ A } there is no undershoot).
Remark 5 It is interesting to notice that, due to the competition between terms a_{32}and b_{32}, the considered Lac Operon model is not a monotone system according to the definition in[16], where a different model was considered.
MAPK signaling pathway
Mitogenactivated protein (MAP) kinases are proteins that respond to the binding of growth factors to cell surface receptors. The pathway consists of three enzymes, MAP kinase, MAP kinase kinase (MAP2K) and MAP kinase kinase kinase (MAP3K) that are activated in series. By activation or phosphorylation, we mean the addition of a phosphate group to the target protein. Extracellular signals can activate MAP3K, which in turn phosphorylates MAP2K at two different sites; in the last round, MAP2K phosphorylates MAPK at two different sites. The MAP kinase signaling cascade can transduce a variety of growth factor signals, and has been evolutionary conserved from yeast to mammals.
Several experimental studies have highlighted the presence of feedback loops in this pathway, which result in different dynamic properties. This work will focus on a specific positivefeedback topology, where doublyphosphorylated MAPK has an activation effect on MAP3K. Such positive feedback has been extensively studied in the literature, since the biochemical analysis of Huang and Ferrell [37, 38] on the MAPK cascade found in Xenopus oocytes. In this type of cells, Mos (MAP3K) can activate MEK (MAP2K) through phosphorylation of two residues (converting unphosphorylated MEK to monophosphorylated MEKP and then bisphosphorylated MEKPP). Active MEK then phosphorylates p42 (MAPK) at two residues. Active p42 can then promote Mos synthesis, completing the closed positivefeedback loop.
where x is concentration of Mos (MAP3K), y_{1} is the concentration of unphosphorylated MEK (MAP2K), y_{2} is the concentration of phosphorilated MEKP, y_{3} is the concentration of MEKPP, z_{1}, z_{2} and z_{3} are respectively the concentrations of unphosphorylated, phosphorylated and doublyphosphorylated p42 (MAPK). Finally, u is the input to the system.
The term μx_{7} introduces the positive feedback loop and represents a key parameter for the analysis to follow. A preliminary screening of the system immediately highlights the following properties. Function b_{11}(x_{1})x_{1}, functions c_{23}(x_{3}), b_{21}(x_{2}), a_{41}(x_{3}) and b_{44}(x_{4})x_{4}, functions c_{56}(x_{6}), b_{54}(x_{5}), a_{74}(x_{6}) and b_{77}(x_{7})x_{7} are increasingasymptoticallyconstant. Moreover, a_{31}(x_{2}) = b_{21}(x_{2}), c_{34}(x_{4}) = b_{44}(x_{4})x_{4}, b_{31}(x_{3}) = a_{41}(x_{3}), b_{33}(x_{3})x_{3} = c_{23}(x_{3}) and a_{64}(x_{5}) = b_{54}(x_{5}), c_{67}(x_{7}) = b_{77}(x_{7})x_{7}, b_{64}(x_{6}) = a_{74}(x_{6}), b_{66}(x_{6})x_{6} = c_{56}(x_{6}). We assume c_{10} to be a positiveconstant.
with k ≐ x_{2}(0) + x_{3}(0) + x_{4}(0) and h ≐ x_{5}(0) + x_{6}(0) + x_{7}(0). Since x_{ i } ≥ 0, all the variables but x_{1} are bounded. The system can be studied by removing variables x_{3} = k  x_{2}  x_{4} and x_{6} = h  x_{5}  x_{7}. We must assume that otherwise no equilibrium is possible. The following result is proved in Section SIV of the Additional File.
Proposition 6 For μ = 0 the system admits a unique globally asymptotically stable equilibrium.
For μ > 0, the system may have multiple equilibria, for specific choices of the involved functions a, b, c.
For μ > 0 suitably large and a_{17}(x_{1}) lower bounded by a positive number, then the system has no equilibria.
For μ > 0 suitably bounded and a_{17}(x_{1}) increasing, or nondecreasing, and bounded, if multiple simple^{2}equilibria exist, then such equilibria are alternatively stable and unstable. In the special case of three equilibria, then the system is bistable.
For μ > 0 suitably bounded and a_{17}(x_{1}) increasing asymptotically unbounded, then the number of equilibria is necessarily even (typically 0 or 2). Moreover, if we assume that there exists μ* > 0 such that the system admits two distinct equilibria for any 0 < μ ≤ μ*, then one is stable, while the other is unstable.
The proof of this last proposition also shows that multiple equilibria x^{ A } , x^{ B } ,... have a partial order: while and have the reverse order and
Remark 6 The simplest case of constant a_{17}has been fully developed in[17]^{3}and[16], and it turns out that the system may exhibit bistability for suitable values of the feedback strength μ. Here we show that, for constant a_{17}, bistability is actually a robust property. Our results are consistent with the fact that the MAPK cascade is a monotone system and some of them could be demonstrated with the same tools used in[16, 17]. With respect to such literature, our contribution is that of inferring properties such as number of equilibria and mono or bistability starting from qualitative assumptions on the dynamics of the model, without invoking monotonicity.
Remark 7 Finally, it is necessary to remark that our results on the MAPK pathway robust behaviors hold true given the model (19) and its structure. Other work in the literature shows that feedback loops are not required to achieve a bistable behavior in the MAPK cascade[38], when the dual phosphorylation and dephosphorylation cycles are nonprocessive (i.e. sites can be phosphorylated/dephosphorylation independently) and distributed (i.e. the enzyme responsible for phosphorylation/dephosphorylation is competitively used in the two steps).
Conclusions
A property is structurally robust if it is satisfied by a class of systems of a given structure, regardless the choice of specific expressions adopted and of the parameter values in the model. We have considered five relevant biological examples and proposed to capture their dynamics with parameterfree, qualitative models. We have shown that specific robust properties of such models can be assessed by means of solid theoretical tools based on Lyapunov methods, setinvariance theory and matrix theory. Robustness is often tested through simulations, at the price of exhaustive campaigns of numerical trials and, more importantly, with no theoretical guarantee of robustness. We are far from claiming that numerical simulation is useless. It it important, for instance, to falsify "robustness conjectures" by finding suitable numerical counterexamples. Furthermore, for very complex systems in which analytic tools can fail, simulation appears be the last resort. Indeed a limit of the considered theoretical investigation is that its systematic application to more complex cases is challenging. However, the set of techniques we employed can be successfully used to study a large class of simple systems, and are in general suitable for the analytical investigation of structural robustness of biological networks, complementary to simulations and experiments.
Notes
^{1}The concentration mismatch is more "softly" controlled, since the derivative of the difference is not influenced by the nonlinear term b_{12}(x_{1}, x_{2})x_{2} = b_{21}(x_{1}, x_{2})x_{1}.
^{2}I.e. the nullclines have no common tangent lines.
^{3}Cf. the erratum: http://www.math.rutgers.edu/~sontag/FTPDIR/angeliferrellsontagpnas04errata. txt and [42].
Declarations
Acknowledgements
The authors acknowledge financial support by the National Science Foundation (NSF) grant CCF0832824 (The Molecular Programming Project). We are grateful to R. M. Murray, for helpful advise and discussions, and to the Reviewers for their constructive comments.
Authors’ Affiliations
References
 Kitano H: Systems biology: A brief overview. Science. 2002, 295 (5560): 16621664. 10.1126/science.1069492View ArticlePubMedGoogle Scholar
 Barkai N, Leibler S: Robustness in simple biochemical networks. Nature. 1997, 387 (6636): 913917. 10.1038/43199View ArticlePubMedGoogle Scholar
 Alon U, Surette MG, Barkai N, Leibler S: Robustness in bacterial chemotaxis. Nature. 1999, 397 (6715): 168171. 10.1038/16483View ArticlePubMedGoogle Scholar
 Yi TM, Huang Y, Simon MI, Doyle J: Robust perfect adaptation in bacterial chemotaxis through integral feedback control. Proceedings Of The National Academy Of Sciences Of The United States Of America. 2000, 97 (9): 46494653. 10.1073/pnas.97.9.4649PubMed CentralView ArticlePubMedGoogle Scholar
 Ma W, Trusina A, ElSamad H, Lim WA, Tang C: Defining Network Topologies that Can Achieve Biochemical Adaptation. Cell. 2009, 138 (4): 760773. 10.1016/j.cell.2009.06.013PubMed CentralView ArticlePubMedGoogle Scholar
 Prill RJ, Iglesias PA, Levchenko A: Dynamic Properties of Network Motifs Contribute to Biological Network Organization. PLoS Biology. 2005, 3 (11): e343.PubMed CentralView ArticlePubMedGoogle Scholar
 Kwon YK, Cho KH: Quantitative analysis of robustness and fragility in biological networks based on feedback dynamics. Bioinformatics. 2008, 24 (7): 987994. 10.1093/bioinformatics/btn060View ArticlePubMedGoogle Scholar
 GómezGardenes J, Y M, Floría LM: On the robustness of complex heterogeneous gene expression networks. Biophysical Chemistry. 2005, 115: 225229. 10.1016/j.bpc.2004.12.033View ArticlePubMedGoogle Scholar
 Gorban A, Radulescu O: Dynamical robustness of biological networks with hierarchical distribution of time scales. IET Systems Biology. 2007, 1 (4): 238246. 10.1049/ietsyb:20060083View ArticlePubMedGoogle Scholar
 Kartal O, Ebenhöh O: Ground State Robustness as an Evolutionary Design Principle in Signaling Networks. PLoS ONE. 2009, 4 (12): e8001.PubMed CentralView ArticlePubMedGoogle Scholar
 Aldana M, Cluzel P: A natural class of robust networks. Proceedings of the National Academy of Sciences of the United States of America. 2003, 100 (15): 87108714. 10.1073/pnas.1536783100PubMed CentralView ArticlePubMedGoogle Scholar
 Tian T: Robustness of mathematical models for biological systems. ANZIAM J. 2004, 45: C565C577.Google Scholar
 Shinar G, Milo R, Rodrìguez Martìnez M, Alon U: Inputoutput robustness in simple bacterial signaling systems. Proceedings of the National Academy of Sciences. 2007, 104: 19931199935. 10.1073/pnas.0706792104.View ArticleGoogle Scholar
 Shinar G, Feinberg M: Structural Sources of Robustness in Biochemical Reaction Networks. Science. 2010, 327 (5971): 13891391. 10.1126/science.1183372View ArticlePubMedGoogle Scholar
 Feinberg M: Chemical reaction network structure and the stability of complex isothermal reactors  I. The deficiency zero and deficiency one theorems. Chemical Engineering Science. 1987, 42: 22292268. 10.1016/00092509(87)800994.View ArticleGoogle Scholar
 Sontag E: Monotone and nearmonotone biochemical networks. Systems and Synthetic Biology. 2007, 1: 5987. 10.1007/s1169300790059PubMed CentralView ArticlePubMedGoogle Scholar
 Angeli D, Ferrell JE, Sontag ED: Detection of multistability, bifurcations, and hysteresis in a large class of biological positivefeedback systems. Proceedings of the National Academy of Sciences of the United States of America. 2004, 101 (7): 18221827. 10.1073/pnas.0308265100PubMed CentralView ArticlePubMedGoogle Scholar
 Jacquez J, Simon C: Qualitative Theory of Compartmental Systems. SIAM Rev. 1993, 35: 4379. 10.1137/1035003.View ArticleGoogle Scholar
 Abate A, Tiwari A, Sastry S: Box Invariance for biologicallyinspired dynamical systems. 46th IEEE Conference on Decision and Control, New Orleans, LA. 2007, 51625167.Google Scholar
 ElSamad H, Prajna S, Papachristodoulou A, Doyle J, Khammash M: Advanced Methods and Algorithms for Biological Networks Analysis. Proceedings of the IEEE. 2006, 94 (4): 832853. 10.1109/JPROC.2006.871776.View ArticleGoogle Scholar
 Radde N, Bar N, Banaji M: Graphical methods for analysing feedback in biological networks  A survey. Int J Syst Sci. 2010, 41: 3546. 10.1080/00207720903151326.View ArticleGoogle Scholar
 Rouche N, Habets P, Laloy M: Stability theory by Liapunov's direct method. 1977, 22: [Applied Mathematical Sciences, xii+396 pp. ISBN 0387902589], New York: SpringerVerlagView ArticleGoogle Scholar
 Blanchini F: Set invariance in control  a survey. Automatica. 1999, 35 (11): 17471767. 10.1016/S00051098(99)001132.View ArticleGoogle Scholar
 Blanchini F, Miani S: Settheoretic methods in control. 2008, 22: [Systems & Control: Foundations & Applications], Boston: BirkhäuserGoogle Scholar
 De Jong H: Modeling and simulation of genetic regulatory systems: a literature review. Journal of Computational Biology. 2002, 9: 67103. 10.1089/10665270252833208View ArticlePubMedGoogle Scholar
 Alon U: An Introduction to Systems Biology: Design Principles of Biological Circuits. 2006, Chapman & Hall/CRCGoogle Scholar
 Nikolov S, Yankulova E, Wolkenhauer O, Petrov V: Principal difference between stability and structural stability (robustness) as used in systems biology. Nonlinear Dynamics Psychol Life Sci. 2007, 11 (4): 413433.PubMedGoogle Scholar
 Mangan S, Zaslaver A, Alon U: The coherent feedforward loop serves as a signsensitive delay element in transcription networks. Journal of Molecular Biology. 2003, 334: 197204. 10.1016/j.jmb.2003.09.049View ArticlePubMedGoogle Scholar
 Levine E, Zhang Z, Kuhlman T, Hwa T: Quantitative Characteristics of Gene Regulation by Small RNA. PLoS Biology. 2007, 5 (9): e229.PubMed CentralView ArticlePubMedGoogle Scholar
 Mitarai N, J B, S K, S S, Z C, E M, Sneppen K: Dynamic features of gene expression control by small regulatory RNAs. Proceedings of the National Academy of Sciences of the United States of America. 2009, 106 (26): 1065510659. 10.1073/pnas.0901466106PubMed CentralView ArticlePubMedGoogle Scholar
 Mehta P, Goyal S, Wingreen NS: A quantitative comparison of sRNAbased and proteinbased gene regulation. Mol Syst Biol. 2008, 4:Google Scholar
 Williamson T, Schwartz JM, Kell DB, Stateva L: Deterministic mathematical models of the cAMP pathway in Saccharomyces cerevisiae. BMC Systems Biology. 2009, 3:Google Scholar
 Jacob F, Perrin D, Sanchéz C, Monod : L'opéeron: groupe de gènes à expression coordonnée par un opérateur. J C R Acad Sci. 1960, 250: 17271729.Google Scholar
 Vilar JMG, Guet C, Leibler S: Modeling network dynamics: the lac operon, a case study. Journal of Cell Biology. 2003, 161 (3): 471476. 10.1083/jcb.200301125PubMed CentralView ArticlePubMedGoogle Scholar
 Yildirim N, Mackey M: Feedback Regulation in the Lactose Operon: A Mathematical Modeling Study and Comparison with Experimental Data. Biophysical Journal. 2003, 84 (5): 28412851. 10.1016/S00063495(03)700137PubMed CentralView ArticlePubMedGoogle Scholar
 Yildirim N, Santillan M, Horike D, Mackey M: Dynamics and bistability in a reduced model of the lac operon. Chaos. 2004, 14 (2): 279292. 10.1063/1.1689451View ArticlePubMedGoogle Scholar
 Huang CYF, Ferrell JJ: Ultrasensitivity in the mitogenactivated protein kinase cascade. Proceedings Of The National Academy Of Sciences Of The United States Of America. 1996, 93: 1007810083. 10.1073/pnas.93.19.10078PubMed CentralView ArticlePubMedGoogle Scholar
 Ferrell J, James E, Machleder EM: The Biochemical Basis of an AllorNone Cell Fate Switch in Xenopus Oocytes. Science. 1998, 280 (5365): 895898. 10.1126/science.280.5365.895View ArticlePubMedGoogle Scholar
 Markevich NI, Hoek JB, Kholodenko BN: Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J Cell Biol. 2004, 164 (3): 353359. 10.1083/jcb.200308060PubMed CentralView ArticlePubMedGoogle Scholar
 Qiao L, Nachbar RB, Kevrekidis IG, Shvartsman SY: Bistability and Oscillations in the HuangFerrell Model of MAPK Signaling. PLoS Comput Biol. 2007, 3 (9): e184.PubMed CentralView ArticleGoogle Scholar
 Ferrell JE, Pomerening JR, Kim SY, Trunnell NB, Xiong W, Huang CYF, Machleder EM: Simple, realistic models of complex biological processes: Positive feedback and bistability in a cell fate switch and a cell cycle oscillator. FEBS letters. 2009, 583 (24): 39994005. 10.1016/j.febslet.2009.10.068View ArticlePubMedGoogle Scholar
 Russo C, Giuraniuc C, Blossey R, Bodart JF: On the equilibria of the MAPK cascade: Cooperativity, modularity and bistability. Physica A: Statistical Mechanics and its Applications. 2009, 388 (24): 50705080. 10.1016/j.physa.2009.08.018.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.