In this section we will analyze five biological networks as case studies. Three of such examples, the L-arabinose, 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 well-known 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 multi-stable 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 L-arabinose network
The arabinose network is analyzed in [28] as an example of feedforward loop. Two genes araBAD and araFGH are regulated by two transcription factors, AraC and CRP. AraC is a repressor, but turns into an activator when bound to the sugar L-arabinose. CRP is an activator when bound to the inducer cyclic AMP (cAMP), which is produced when cells are starving upon glucose (not produced during growth on glucose). CRP also binds to the araC promoter and enhances transcription of AraC, which has a significant basal rate of expression (i.e. it is produced by the cell also in absence of inducer CRP). A very simple model for this network can be derived by defining the state variables x1 and x2, respectively the concentrations of the transcription factor AraC and of the output protein araBAD. The concentration of the transcription factor CRP is considered an external input u:
where α1, α2 are the degradation and dilution rates of x1, x2 respectively. The basal production rate of x1 (AraC) is p1. The activation pathways are modeled by Hill functions f (u, K) = uH /(KH + uH ), where H is the Hill coefficient and K
ij
are the activation thresholds. The model can be recast into the general structure (4), which includes model (5) as special case::
where u is nonnegative-constant, c1, b11 and b22 are positive-constant, while c1u(u) and c2u 1(u) are sigmoidal with respect to u, the latter increasing with respect to x1. The graph representation of this network is in Figure 3A.
For this elementary network the analysis is straightforward. Variable x1 is not affected by x2. Since c1u(u) is bounded, x1 is also bounded and converges to an equilibrium point
which is monotonically increasing in u. In turn, x2 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 c1 prevents x1(t) and x2(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
Small regulatory RNAs (sRNA) play a fundamental role in the stress response of many bacteria and eukaryotes. In short, when the organism is subject to a stimulus that threatens the cell survival, certain sRNA species are transcribed and can down-regulate the expression of several other genes. For example, when E. coli cells are lacking a source of iron, the sRNA RyhB (normally repressed by the ferric uptake regulator Fur) is expressed and rapidly induces the degradation of at least other 18 RNA species encoding for non-essential proteins which use up Fe molecules. This allows essential iron-dependent pathways to use the low amount of available iron. Quantitative studies of the sRNA pathways have been carried out in [29–31]. Let us define x1 as the RNA concentration of the species which is targeted by the sRNA and x2 as the concentration of sRNA. The model often proposed in the literature is:
where α1, α2 are the transcription rates of x1 and x2 respectively, β1, β2 are their degradation rates (turnover), and k is the binding rate of the species x1 and x2. The formation of the inactive complex x1 · x2 corresponds to a depletion of both free molecules of x1 and x2. If α1 < α2 the pathway successfully suppresses the expression of the non-essential gene encoded by x1. This model can be embedded in the general family:
by setting b12 = kx1 and b21 = kx2 (note that b12(0) = b21(0) as required). From our list of properties: c1, c2, b11 and b22 are positive-constant; b12(x1, x2) and b21(x1, x2) are increasing-asymptotically-unbounded in both variables; and b12(x1, x2)x2 = b21(x1, x2)x1 at all times. This network can be represented with the graph in Figure 2A.
The sRNA system is positive, because the nonnegativity Assumptions 1 and 4 are satisfied. The preliminary screening of this network tells us that each variable produce an inhibition control on the other, which increases with the variable itself. In other words x1 is "less tolerant" to an increase of x2 if the latter is present in a large amount. This means that the sum x1 + x2 is strongly kept under control. Also the mismatch between the two variables is controlled. 1 To prove stability of the (unique) equilibrium
, we will use the 1-norm as Lyapunov function
(see Figure 2B). This choice has a remarkable interpretation: denoting by
and
the sum and the difference of the two variables (referred to the equilibrium) we have
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.
Proposition 2 The variables of system (8) are bounded for any initial condition x1(0), x2(0) ≥ 0. The system admits a unique asymptotically stable equilibrium point
and the convergence is exponential:
for some β > 0 and any x1(0) ≥ 0, x2(0) ≥ 0. Moreover, no oscillations are possible around the equilibrium, in the sense that the condition
or
occurs at most once.
Proof To prove boundedness of the variables we need to show the existence of an invariant set
Proposition 1 guarantees that the positivity constraints are respected. Then we just need to show that the constraint x1 + x2 ≤ κ cannot be violated for sufficiently large κ > 0. The derivative of function s (x1, x2) = x1 + x2 is
Define κ = (c1 + c2)/min {b11, b22} then for s(x1, x2) > κ the derivative becomes negative so s(x1, x2) cannot exceed κ (See Theorem 1).
Boundedness of the solution inside a compact set assures the existence of an equilibrium point. Let
be any point in which the following equilibrium conditions holds:
The behavior of the candidate Lyapunov function:
will be examined in the different sectors represented in Figure 2B. Let us start by considering the sector defined by
and
(APB in Figure 2B) for which
. In such a sector the Lyapunov derivative is:
where we have subtracted the null terms (10) and where we have exploited the fact that b12(x1, x2)x1 = b21(x1, x2)x2 is increasing in both variables. The inequality (CPD in Figure 2B)
can be similarly proved to hold in the sector
and
.
Consider the sector defined by
and
(DPA in Figure 2B) for which
in such a sector the Lyapunov derivative is:
Note that in the last step we have added and subtracted the null terms (10). In the opposite sector (BPC in Figure 2B)
and
, we can prove that
. We just proved that
with β = min{b11, b22}. 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.
We can apply Nagumo's theorem: consider the half-line PA originating in P, where
and
. Therefore we have that (again by adding the null term in (10)):
Similarly, on half-line PD where
and
, by considering (10) we derive
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 c1(t) and c2(t) are time-varying, but bounded, we have convergence to a neighborhood whose amplitude, obviously, depends on the bounds of c1(t) and c2(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 up-regulate entire groups of transcripts and therefore increase or decrease c2(t).
The following corollary evidences the positive influence of c2, which is positive over x2 and negative over x1.
Corollary 1 Assume that x1(0), x2(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.
Proof The steady state values
and
are respectively monotonically decreasing and increasing functions of c2. Indeed, consider the steady-state condition
and regard it as a differentiable map (x1, x2) → (c1, c2). By the uniqueness proved in Proposition 2 the map is invertible. The Jacobian of the inverse map is the inverse of the Jacobian
where
(keep in mind that b21(x1, x2)x1 = b12(x1, x2)x2). The sign of the entries in the second column are negative and positive respectively, therefore, the steady-state values
and
are decreasing and increasing functions of c2.
The absence of overshoot and undershoot is an immediate consequence of the invariance of the sector
and
previously proved.
Obviously, decreasing c2 increases x1 and decreases x2 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
The cyclic adenosine monophosphate (cAMP) pathway can activate enzymes and regulate gene expression based on sensing of extracellular signals. Such signals are sensed by the G protein-coupled receptors on the cell membrane. When a receptor is activated by its extracellular ligand, a series of conformational changes are induced in the receptor and in the attached intracellular G protein complex; the latter activates adenylyl cyclase, which catalyzes the conversion of ATP in cAMP. In yeast, cAMP causes the activation of the protein kinase A (PKA), which in turn regulates the cell growth, metabolism and stress response. Stochastic models are usually proposed for numerical analysis of the cAMP pathway. However, the cAMP pathway components in yeast are present in high numbers and a deterministic modeling approach is adopted in [31]. In such paper, the pathway is decomposed in several modules, here we consider the simplified cAMP Model A, which focuses on the parts of the pathway best characterized in the literature:
where x1 is the concentration of active G protein, x2 is the concentration of active PKA protein, x3 is the concentration of cAMP and u is the concentration of glucose input to the network. The parameters
and
model the "feedback" effect introduced by two phosphodiesterases (Pde1p and Pde2p) on the cAMP concentration. The numerical analysis in [31] typically shows that the cAMP concentration (x3) responds with a large overshoot to steps in the glucose (u, input) concentration. We will analytically explore the dynamic behavior of x3, showing that under certain assumptions, a bounded overshoot is a robust characteristic in the system. The parameters k
F
and k
R
are forward and reverse reaction rates for the formation of active x1 and x2. Mass conservation allows to express the active protein amounts as a function of the total number of molecules,
. The nonlinear expressions in equation x3 are derived by Michaelis-Menten enzymatic steps. We can re-write the above equations according to the general model (4):
where u is the external signal and where b23 = 0 for x2 = 0 and b32 = 0 for x3 = 0. A qualitative graphical representation of this network is in Figure 3B.
Our preliminary analysis allows us to assume: a1u, a23: decreasing-exactly-null with threshold values
and
; d32, a31: decreasing-asymptotically-null, b32 and g33 = b33(x3)x3: increasing-asymptotically-constant; b11, b22are positive-constant.
It is immediate to notice that for constant u, x1 robustly converges asymptotically to its equilibrium value such that
The solution
of the previous equation is uniquely defined for each u since the function ξ-1(x1) on the right is strictly increasing and grows to infinity, precisely
. Biologically, this means that external glucose signals are mapped to internal active G-protein concentration with a bijection, before saturating.
Also we have to note that the model is consistent with mass conservation: since a1u(x1) and a23(x2) are zero above the thresholds
and
, we have that
and
for
and
, respectively; therefore we assume
,
, for all t ≥ 0.
Proposition 3 There exists an equilibrium for system (12) if and only if
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 x1 is at its equilibrium value
. Furthermore, under a suitable condition a performance bound on the transient values of x3(t) can be given.
Proposition 4 Assume that x1has reached its steady state
. Then, the unique equilibrium point is globally attractive for any initial condition x2(0), x3(0) ≥ 0. Moreover, assume that
then we can give the following bound for the transient of x3(t)
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 b33(x3)x3at "full force" (x3suitably large) dominates the activating term d32(x2) + a31(x2)ξ when x2is 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 x3(0), may exhibit a spike of cAMP x3which is bounded by (15), if condition (14) is satisfied. If x
3
(0) is small, then the bound is d32(0) + a31(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 x3(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 multi-stable.
The state variables of the ODE model we will study are the concentration of nonfunctional permease protein x1; the concentration of functional permease protein x2; the concentration of inducer (allolactose) inside the cell x3, and the concentration of β-galactosidase x4, a quantity that can be experimentally measured. The concentration of inducer external to the cell is here denoted as an input function u.
where β1, β2, δ1, δ2, δ3 and γ are constants and f
i
are functions that are experimentally measurable. In particular, at low inducer concentrations,
, where k
i
's are constant; at high x3 concentrations f1 saturates. The functions f2 and f3 are assumed to depend hyperbolically on their arguments. According to the proposed setup, the previous equations can rewritten as follows:
where c13(x3) = f1(x3), b11 = δ1, a21 = β1, b22 = δ2, a32(u) = f2(u) =, b32(x3) = f3(x3), c3u= β2, b33 = δ3, c43(x3) = γ f1(x3) and b44 = δ4. This corresponds to the network in Figure 3C.
From our preliminary analysis step: c13 is constant-sigmoidal, a32(u) and b32(x3) are increasing-asymptotically-constant, and the remaining functions a21, b11, b22 and b33 are positive-constant.
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 x4 does not affect any other chemical species: therefore, the fourth equation can be considered separately. As long as the inducer concentration of x3 within the cell reaches an equilibrium
, x4 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 S-III of the Additional File), neglecting the linear term c3uu 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 x1 receives a bounded signal from x3 and the degradation term -b11x1 keeps x1 bounded. In turn, x2 remains bounded. The inducer concentration x3 receives a bounded signal form u and x2; therefore x3 stays bounded as well, being both a32(u) and b32(x3) 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 Hill-type 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 xA , xB , xC ,... ∈ IR3(typically three) for intermediate values of u. If multiple equilibria exist, then they are ordered in the sense that xA ≤ xB ≤ xC ... 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, xA , xB , xC 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 S-III of the Additional File. The cone-invariance property implies that the state variables cannot exhibit oscillations around their equilibria. For instance, if xA is the first (hence stable) equilibrium, given any initial condition upper bounded by xA (x(0) xA ) in the domain of attraction, the convergence to xA has no overshoot (and if x(0) ≥ xA there is no undershoot).
Remark 5 It is interesting to notice that, due to the competition between terms a32and b32, 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
Mitogen-activated 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 positive-feedback topology, where doubly-phosphorylated 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 MEK-P and then bisphosphorylated MEK-PP). Active MEK then phosphorylates p42 (MAPK) at two residues. Active p42 can then promote Mos synthesis, completing the closed positive-feedback loop.
The presence of such positive feedback in the MAPK cascade has been linked to a bistable behavior: the switch between two stable equilibria in Xenopus oocytes denotes the transition between the immature and mature state. A standard ODE model for the cascade is proposed in [17], where the authors demonstrate bi-stability of the system by applying the general theory of monotone systems. We adopt such model, which is reported below:
where x is concentration of Mos (MAP3K), y1 is the concentration of unphosphorylated MEK (MAP2K), y2 is the concentration of phosphorilated MEK-P, y3 is the concentration of MEK-PP, z1, z2 and z3 are respectively the concentrations of unphosphorylated, phosphorylated and doubly-phosphorylated p42 (MAPK). Finally, u is the input to the system.
While bi-stability may occur due to other phenomena, such as multisite phosphorylation [39], rather than due to feedback loops, a large body of literature focuses on bi-stability induced by the positive-feedback in the Huang-Ferrel model in Xenopus[40, 41] reported above. In [37] the feedback f (u) was characterized, through in vitro studies, as an activating Hill-function with high cooperativity. In [17] instead, f (u) was assumed to be a first order linear term in the concentration of active MAP3K, x7. In Proposition 6, we will explore the effects of different qualitative functional assumptions on the feedback loop dynamics f (u). We will highlight that the system loses its well-known bi-stability not only in the absence of feedback, but also when the feedback becomes unbounded. An unbounded positive feedback would be caused, for instance, by an autocatalytic process of MAP3K activation, mediated by active MAPK. We choose to rewrite the above model as follows:
The term μx7 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 b11(x1)x1, functions c23(x3), b21(x2), a41(x3) and b44(x4)x4, functions c56(x6), b54(x5), a74(x6) and b77(x7)x7 are increasing-asymptotically-constant. Moreover, a31(x2) = b21(x2), c34(x4) = b44(x4)x4, b31(x3) = a41(x3), b33(x3)x3 = c23(x3) and a64(x5) = b54(x5), c67(x7) = b77(x7)x7, b64(x6) = a74(x6), b66(x6)x6 = c56(x6). We assume c10 to be a positive-constant.
The graph in Figure 3D can be partitioned considering three aggregates of variables, precisely {x1}, Σ234 = (x2, x3, x4) and Σ 567 = {x5, x6, x7}. Signal x1 is the only input for Σ234, signal x4 is the only input for Σ567. Then x7 is fed back to the first subsystems by arc a17. Without the positive feedback loop, we will demonstrate that the system is a pure stable cascade. Note also that Σ234 and Σ567 can be reduced since
and
and therefore the following sums are constant
with k ≐ x2(0) + x3(0) + x4(0) and h ≐ x5(0) + x6(0) + x7(0). Since x
i
≥ 0, all the variables but x1 are bounded. The system can be studied by removing variables x3 = k - x2 - x4 and x6 = h - x5 - x7. We must assume that
otherwise no equilibrium is possible. The following result is proved in Section S-IV 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 a17(x1) lower bounded by a positive number, then the system has no equilibria.
For μ > 0 suitably bounded and a17(x1) increasing, or non-decreasing, and bounded, if multiple simple2equilibria 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 a17(x1) 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 xA , xB ,... have a partial order:
while
and
have the reverse order
and 
Remark 6 The simplest case of constant a17has been fully developed in[17]3and[16], and it turns out that the system may exhibit bi-stability for suitable values of the feedback strength μ. Here we show that, for constant a17, bi-stability 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 bi-stability 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 de-phosphorylation cycles are non-processive (i.e. sites can be phosphorylated/de-phosphorylation independently) and distributed (i.e. the enzyme responsible for phosphorylation/de-phosphorylation is competitively used in the two steps).