Dynamics and Structure in Cell Signaling Networks: Off-State Stability and Dynamically Positive Cycles

The signaling system is a fundamental part of the cell, as it regulates essential functions including growth, differentiation, protein synthesis, and apoptosis. A malfunction in this subsystem can disrupt the cell significantly, and is believed to be involved in certain diseases, with cancer being a very important example. While the information available about intracellular signaling networks is constantly growing, and the network topology is actively being analyzed, the modeling of the dynamics of such a system faces difficulties due to the vast number of parameters, which can prove hard to estimate correctly. As the functioning of the signaling system depends on the parameters in a complex way, being able to make general statements based solely on the network topology could be especially appealing. We study a general kinetic model of the signaling system, giving results for the asymptotic behavior of the system in the case of a network with only activatory interactions. We also investigate the possible generalization of our results for the case of a more general model including inhibitory interactions too. We find that feedback cycles made up entirely of activatory interactions (which we call dynamically positive) are especially important, as their properties determine whether the system has a stable signal-off state, which is desirable in many situations to avoid autoactivation due to a noisy environment. To test our results, we investigate the network topology in the Signalink database, and find that the human signaling network indeed has only significantly few dynamically positive cycles, which agrees well with our theoretical arguments.


Introduction
Biochemical interactions in a cell form complex networks with its elements carrying out many different tasks [1]. An important part is the signal transduction system which allows the cell to adapt to the environment. The working of this signaling system has been actively studied, both theoretically and experimentally, with databases now being available with growing information about the specific proteins making up the signal transduction network and their interactions [2][3][4][5][6]. As the signaling system is an essential part of the cell, a disruption in its functioning is believed to be a significant factor in many diseases, especially in cancer [7][8][9]. Thus, a better understanding of the signaling system will possibly allow more effective treatments. Also, from a theoretical point of view, the study of cell signaling can reveal how biochemical interactions are able to form a complex self-organizing system, capable of processing information in a noisy environment.
The signaling system is naturally modeled as a directed graph: the nodes are the protein species possibly present in a cell, and an edge points from node i to node j if protein species i affects protein species j. A kinetic model is usually a set of ordinary differential equations where the independent variables represent the concentration of the possible states of the proteins. An edge in the graph implies that the some of the variables associated with protein i is present on the right-hand-side of the equations concerning some of the variables associated with protein j.
In this paper we investigate the general mathematical model introduced by Heinrich et al. [10,11]. We follow the simple hypothesis that activation in the signaling system should only occur as a result of an external or internal signal (e.g. a receptor is activated) [11,12]. While in some cases this is not true (a signal only tunes the output [13]), this might prove as a correct basic principle for the part of the system which has to respond to the stimuli. Mathematically, this means that the signal-off state of the dynamical system describing the signaling network has to be stable and also attractive in the absence of stimulation so that the system will relax to it, and will not be spontaneously activated by any noise present [12].
The original model of Heinrich et al. includes only positive interactions: a protein species activating other proteins. In this case, we identify the global attractor of the system, which is reached in the presence of any constant external inputs. The addition of inhibitory interactions gives rise to more complex behavior generally, but the property of whether the network has a stable signal-off state is preserved. To assess our theoretical arguments, we investigate the human signaling network found in the Signalink database [4], and find that having a stable signal-off state indeed seems to be an important factor shaping the network topology.

A Network Model with Positive Links
A possible model of the signaling network is that of Heinrich et al. [10,11]. In this model, each protein has two possible states (active and inactive, activation is achieved via phosphorylation), with only proteins in their active state interacting with others. A protein can become activated by another protein or by a stimulated receptor which is considered to be the external input driving the system. Deactivation occurs via phosphatases, which are considered to work independently of the rest of the system. This model is general in the sense that multiple interconnected signaling pathways can be included in one network, while it contains only simple activating interactions. This means that this model is targeted at the propagation of some signal (primarily considering, but not restricted to an external stimulus getting into the cell's nucleus) not including some advanced features like adaptation which require a more specialized treatment [6,13].

The Equations Considered
Using mass action kinetics, the system can be described by the following ordinary differential equations: wherex x i and e x x i are the concentrations of the i-th protein in the active and inactive state respectively. The b i constants describe the deactivation of the i-th protein with its rate being independent of the other protein concentrations. The elements of the matrix e a a ij describe the activation of protein i by protein j. We assume that all b i are positive, all e a a ij are nonnegative and the diagonal elements are all zero: e a a ii~0 . The term S i (t) is an external signal affecting the i-th protein, with S i being positive for proteins receiving input, and zero for all others. We will assume that the external signal is constant or changes on a time scale significantly larger than the rest of the system, and write S i (t)~s i . We also make the assumption that the total concentration of each protein species is constant:x x i ze x x i~Ci , which allows us to write Eq. (1) in a dimensionless form: where we defined x i :x x i =C i and a ij :C j e a a ij . Note that now all x i [½0,1. Also, the new matrix elements a ij are scaled by the total protein concentration C j , meaning that the change of total concentrations (e.g. by the synthesis of proteins) can be represented in this model by varying the elements of the matrix a ij . We chose not to include the effects of concentration changes in the model itself, but this allows the comparison of two systems where the C j concentrations differ. In the following, we will denote the vector with components x i by x. Throughout this paper we interpret relations between vectors component-wise, e.g. x §y (or xwy) means that x i §y i (x i wy i respectively) for all i.
If we have no external signal (all s i are zero), the off-state (when x i :0 for all i) is always a stationary solution. In the special case, when the network contains no cycles (i.e. a cascade), the off-state is always stable and globally attractive; a finite-time signal has a response which relaxes to the off-state exponentially [10]. If the network contains cycles, the off-state can be unstable and a nonzero solution of Eq. (2) can persist even in the absence of an external signal. Linear stability of the off-state with respect to the network topology was investigated by Kartal and Ebenhöh [12]. We now give the general attractor of this system in the presence of arbitrary constant inputs.
The relation of Eq. (2) with the interaction graph of the protein network is quite natural. An edge i?j implies that x i is present in the equation concerning _ x x j . The a ij matrix in Eq. (2) can be considered as the appropriately weighted adjacency matrix of the interaction graph.

The General Attractor
We can make a general statement about the solutions of Eq. (2) which hold for any possible combination of coefficients a ij and b i with any combination of stationary inputs s i . In this section we summarize our results. Proof of our statements is presented in the Methods section.
Eq. (2) always has exactly one stable stationary solution which is globally attractive with respect to the valid range of initial conditions x i [½0,1, and may have other stationary solutions which are unstable. Despite that Eq. (2) describes a nonlinear system, no more complex behavior is possible; the system started from any initial condition which is not itself an unstable stationary solution relaxes to the unique stable stationary solution.
The main idea is that if there is a linearly stable stationary solution of Eq. (2), then a Lyapunov function can be used to show that it is attractive for any initial condition from the x §x Ã region. This implies that there can be no other stable stationary solutions as their basins of attraction would overlap. Considering the region where for at least one coordinate x i vx Ã i , we use the properties of the Jacobian matrix J ij~Li _ x x j of the system to prove that regardless of the initial condition (except for a possible unstable stationary solution) the system will reach the x §x Ã region where the Lyapunov function is applicable. Also, using the properties of the Jacobian matrix, we prove the existence of at least one linearly stable stationary solution.
Our results extend the those of Kartal and Ebenhöh [12]; considering the system described by Eq. (2), they examined linear stability and concluded that a linearly stable off-state is a significant property of real-world signaling networks. Using our results, it follows that a stable stationary solution is globally attractive; if the off state is linearly stable, the system will eventually return to it after any external signal has been shut off.

Networks with Inhibitory Interactions
We examined the generalization of Eq. (2) to include inhibitory interactions. We considered the case when a protein in its activated state can deactivate other proteins, which can be described by the following equations: Here the elements of the c ij matrix describe the deactivation of protein i by protein j. We assume that all c ij are nonnegative with the diagonal elements being zero. If we set all c ij~0 , we get back the system (2). Computing the Jacobi matrix for the off-state (considering the case when the external signal s i is zero), we get that it does not include the c ij matrix. This means that the linear stability of the off-state depends only on the elements of the matrix a ij and the deactivation rate constants b i . If the off-state for some particular system with only positive interactions is stable, then it remains stable if we add inhibitory interactions. This is also true for the global attractive properties of the offstate; stability of the off-state implies that it is globally attractive even in the presence inhibitory interactions. The Lyapunovfunction defined for the system with only positive interactions can be applied in this case too. Computing the derivative of this Lyapunov-function, we get that the terms containing the inhibitory interactions are nonpositive for all possible x, meaning that if a Lyapunov-function was valid in a network with positive interactions, then it remains valid after the addition of negative interactions (see the Methods section, especially Eqs. (8) and (15)).
A consequence of this property is that if we are interested in the stability of the off-state of a particular network, then we have to consider only the positive interactions between its nodes; we can restrict ourselves to the analysis of the subgraph which includes the positive links only. Behavior of the system in the absence of an external signal depends on the properties of the strongly connected components of that subgraph.
In the case when we have an external signal or the off-state is unstable, we cannot make a general statement about the solutions of Eq. (3). Numerical simulations show cases of monostability, multistability and periodic solutions, suggesting that the inclusion of inhibitory interactions might be necessary to account for a more complex behavior than which can result from Eq. (2), including a wide range of possible input-output patterns [14].

General Two-state Networks
A more general model of the signaling network would have the following form: Eq. (4) is a possible generalization of Eq. (3). In this case too, every protein has two states, and only proteins in the active state can activate or deactivate other proteins. A i denotes the set of nodes which can activate node i and C i is the set of nodes which can deactivate node i (i.e. there is a positive or negative edge between nodes i and j[A i or C i ). We assume, that for each (i,j) pair, node j belongs to at most one of the sets A i and C i . We also restrict the activation term f i to be nonnegative and a monotonically decreasing function of x i , and monotonically increasing function of all x j [A i . The deactivations term g i is also nonnegative, and a monotonically increasing function of x i and all x j [C i . We also require the h i term describing the spontaneous deactivation of protein i to be positive for all x i =0, thus the concentration of the proteins relaxes to 0 if there are no activating interactions. The term V i represents the external input on protein i. Of course, to ensure that the solutions stay in the (these conditions will arise naturally if we derive Eq. (4) from some kinetic model where the total protein concentrations are preserved). Generally, we do not need the f i , g i , h i , V i functions to be continuous or differentiable; if however, we want to carry out the linear stability analysis of a fixed point, then differentiability is required at least in some neighborhood of the fixed point. The most important consequence of the conditions prescribed on Eq. (4) is that the concentration of each protein relaxes to 0 if it receives no activatory interactions or external input.
Giving the general solutions of a model like Eq. (4) can prove very difficult, with chaotic solutions possible in many cases [15]. This is why methods where a statement about the solutions of the model can be made based only on the network topology can prove very useful.

Consequences of the Network Topology
The structure of Eq. (4) gives a very simple requirement for the system to have a nonzero solution in the absence of an external signal. We can notice that the only positive term on the righthand-side of our equations is that describing activation, i.e. the term corresponding to positive links between proteins. If there are no cycles made up entirely of positive links (in the following sections we shall refer to these as dynamically positive cycles), then the system will relax to its off-state if there is no external signal. This is the generalization of the results we got for the specific model described by Eq. (3). The main assumption for this to be true is that in the absence of activatory interactions, each protein relaxes to its inactive state (e.g. via dephosphorylation by phosphatases). The requirement that a nonzero solution in the absence of an external signal can only persist if there are dynamically positive cycles holds true even if we include more complex processes, e.g. multiple phosphorylation.
On the other hand, a general network will possibly contain dynamically positive cycles. In this case, the stability of the off-state will depend on the nature of the functions involved in Eq. (4). Still, using the property that in the absence of activatory interactions, each protein relaxes to its inactive state, we can conclude that only the strongly connected components (SCCs) formed entirely by positive links (i.e. dynamically positive SCCs) are relevant. This means that instead of examining the properties of the whole network, we only have to evaluate the subnetworks which contain the points in each dynamically positive SCC. For the whole system to have a stable off-state, each of these subnetworks is required to have a stable off-state. This is the extension of the ideas in [12] to networks which include both activatory and inhibitory interactions. As every dynamically positive cycle can possibly give rise to an unstable off-state, having as few such cycles as possible (while retaining the biological functionality of the network) might be desirable to avoid an unwanted activation due to variations in parameter values (e.g. the total concentration of the signaling proteins depends on the expression level of the corresponding genes). In the following section, we test this hypothesis on the human signaling network in the Signalink database.

Analysis of the Signalink Network
The Signalink database contains information about proteins and interactions present in three species: Caenorhabditis elegans, Drosophila melanogaster and Homo sapiens, which were manually collected from the literature [4,16]. The interactions are labeled as either activation or inhibition. The networks of both C. elegans and D. melanogaster show a very simple topology with only a few cycles in them, while the human network has a largest SCC of 81 proteins and a total of 777041 cycles (see table 1). We carried out further analysis on the largest SCC of the human network, measuring the number of dynamically positive cycles. We compared the results with the same computed for random networks generated from the original database. We used three methods for network generation. In method (a) we exchanged the sign among randomly chosen pairs of edges, keeping the start and end nodes on both edges. Note that this method preserves the total number of positive or negative edges (if both edges had the same sign, no alteration was made). In method (b) we exchanged the endpoints of randomly chosen pairs of edges, thus also affecting the network topology (but in this case keeping the sign of the edges). Method (c) is an extension of method (b), where we also exchanged the sign of the chosen pair of edges with probability p~0:5.
Based on our theoretical arguments, we expected the real network to have as few dynamically positive cycles as possible. With respect to the stability of the off-state, having no such cycles would be ideal, while this might not be acceptable as dynamically positive cycles could serve a specific purpose. The network in the Signalink database had 13 dynamically positive cycles. We compared this value with the number of such cycles in 200000 random graphs for each method, where each random graph was generated by making 10000 exchanges in the original network. For all three methods, the distribution of the number of dynamically positive cycles is fat-tailed; number of such cycles spans more than five orders of magnitude. The percentile plot of the measured distribution shows that the majority of the random graphs generated has more dynamically positive cycles than the actual value (see Fig. 1, pay attention to the logarithmic scaling of the y-axis). The ratio of graphs having 13 or less dynamically positive cycles is 1:72% for method (a), 0:039% for method (b) and 1:7525% for method (c). This suggests that the number of dynamically positive cycles in the real network is significantly low. This implies, that the number of such cycles is indeed an important property of the signaling network, and while having some such cycles might not be avoidable, a small number is desirable.
From a theoretical view, the importance of dynamically positive cycles arises from two preliminary hypotheses; the assumption that each protein in the network relaxes to its inactive state if it does not receive activatory interactions; and the requirement that the signaling system has a stable off-state. Finding only a significantly small number of dynamically positive cycles in the human signaling network suggests that these hypotheses are valid for a large part of the system.

Dynamically Positive Cycles in the Human Network
In the case of the human signaling network in the Signalink database, we identified 13 dynamically positive cycles, which form 5 SCCs (three of which consist of only two nodes each). These 5 SCCs are all connected in one way with path made up of positive interactions. Table 2 lists the Uniprot IDs [17] of the proteins involved in them, while Fig. 2 shows the interactions between them.

Metzler Matrices
Let A, a square matrix, with elements A ij , be a Metzler matrix, meaning that A ij §0 for all i=j, while the diagonal elements can be negative [16]. Let us consider the case when at least one element of the diagonal is negative. Let a be the minimum of the diagonal: a: min i A ii v0. Then the matrix A can be written in the form A~b A A{DaDI, where b A A is nonnegative and I is the identity matrix. The matrices b A A and A have the same eigenvectors with the corresponding eigenvalues only being shifted by a. For the matrix b A A we can apply the Perron-Frobenius theorem [17] getting that b A A has a real largest eigenvalue (l l 1 §0 such that Dl l i Dƒl l 1 for all other eigenvalues) with the corresponding eigenvector having only nonnegative elements. For the original matrix A, we have l 1~l l 1 {DaD, and for the real part of all other eigenvalues Rel i ƒl 1 . In the following sections we will refer to l 1 as the largest eigenvalue.
In the case of Eq. (2), the Jacobian matrix is a Metzler matrix: the diagonal elements are negative, while the off-diagonal elements J ij are positive when protein j affects protein i, and zero otherwise. While in a general setting J will not be symmetric, the Perron-Frobenius theorem guarantees that the eigenvalue with the largest real part will be real, and the corresponding u eigenvector can be chosen to have (real) nonnegative elements. In the case, when the signaling network is strongly connected (i.e. every node can be reached via directed links from every other node) the corresponding J matrix will be irreducible. In this case the Perron-Frobenius theorem also states that u will have only positive elements [17].

The Derivative of the Lyapunov Function
Let us consider Eq. (2) and presume that there exists a stationary solution x Ã which is linearly stable. We will prove the existence of at least one such stationary state in the next section. Let J Ã denote the Jacobian matrix computed at the stationary solution: Linear stability means that all eigenvalues of J Ã have negative real parts. A Lyapunov function can be constructed showing that for all initial conditions x (0) §x Ã the system relaxes to the stationary solution x Ã . First let us rewrite Eq. (2) for new variables y:x{x Ã : Notice that if y §0 for some t 0 then y §0 will hold for all t §t 0 . Since the solutions of Eq. (6) are continuous functions of t, the system can only leave the y §0 region if at some time t 1 , at least one component is zero. Substituting y i~0 into Eq. (6) and assuming that y j §0 for all j=i, we will have a nonnegative derivate (recall that all a ij §0 and 0vx Ã i v1): This allows us to define a Lyapunov function which is valid in the y §0 region (note that this will imply that our Lyapunov function is valid for all x §x Ã , not just some neighborhood of it). The form of this Lyapunov function is: where the A i positive constants need yet to be chosen in an appropriate way. We have L(y) §0 and L(y)~0 u y~0. We need to choose the A coefficients in a way which gives dL=dtƒ0 and dL=dt~0 u y~0. Computing the time derivative of L we get: where For all y §0, f (y) is nonnegative. This allows us to disregard f (y) and only focus on the first part of the derivative, which can be written in a simpler way using the Jacobian matrix in the stationary solution defined in Eq. (5): where in parentheses the vector of coefficients A i is multiplied by the transpose of the Jacobian matrix J Ã . We get a negative result for all possible y §0 if all components are negative.
In the case of a strongly connected network, using that J ÃT is a Metzler matrix, we get that the eigenvector corresponding to its largest eigenvalue can be chosen to have only positive elements (see the previous section). If the stationary solution x Ã is linearly stable, the largest eigenvalue of J ÃT is negative. In this case, choosing the A vector of coefficients to be the corresponding eigenvector, we have Aw0 and J ÃT Av0. This choice gives a Lyapunov function which is valid for all y §0.
In the case of a network with multiple SCCs, we can choose the A vector to be a linear combination of the eigenvectors of the submatrices corresponding to each separate strongly connected component. For this, let us reorder the x i independent variables in such way, that the SCCs form blocks in J Ã . In a network with n SCCs, we get a matrix with the following structure: Here J i is the Jacobian matrix of the i-th SCC at the stationary solution x Ã and J ij represents interaction between the i-th and j-th SCC (e.g. the second SCC affects the first if J 12 has nonzero elements). All elements above the diagonal blocks are zero (links between two distinct SCCs can be only in one direction). The stationary solution x Ã requires that all J T i have only eigenvalues with negative real parts. Combining this with the fact that all J T i are now irreducible Metzler-matrices, we get that each submatrix will have a real largest eigenvalue l i v0 and corresponding real eigenvector u i w0. Now we can choose the A vector in Eq. (8) to be a ''linear combination'' of these eigenvectors: A~a 1 u 1 ,a 2 u 2 , . . . a n u n ð Þ . Of course we need to have a i w0 for all i. With this, we have: The first block in the resulting B vector will be negative for any a 1 w0. All other blocks can be separated into two terms with opposite sign: The summation will yield a nonnegative result, and the last term will be negative. Starting from j~2 and heading consecutively to j~n, we can always choose a j constants that are big enough for B j to be negative. In this way, we have constructed a suitable A vector of coefficients for the Lyapunov function in Eq. (8).
Using the properties of the Lyapunov function, we get that for any initial condition x §x Ã the system will converge to its stationary solution x Ã . This also means that there can be at most one stable stationary solution; if there were more, their basins of attraction would overlap.
Note that our arguments presented for the networks with multiple SCCs do not actually require that we know the eigenvalues of J Ã ; knowing that each submatrix J i has only negative eigenvalues is sufficient. Of course, in a general setting the submatrices are not independent: J i will be affected by J j possibly for any jvi. A rather special case is that when we have subsystems that all have stable off-states: any arrangement of these subsystems which does not include feedback loops will also have a stable offstate.

The Lyapunov-function for the Off-state in the Presence of Inhibitory Interactions
Here we consider the special case of Eq. (3) when all s i~0 and the off-state is linearly stable (i.e. all eigenvalues of the Jacobian matrix are negative). We can define the Lyapunov function for this system in a way similar to the previous section: where we require all A i to be positive. Computing the time derivative we get a form similar to Eq. (9): where J Ã is the Jacobian matrix in the off-state, and f is nonnegative, and the inhibitory interactions only appear in f . Using the same argument as in the previous section, we can choose an appropriate Aw0 vector which will give a positive Lyapunov function and a negative time derivative. This means that the system will reach its stationary solution x~0.

The Existence of a Stable Stationary Solution and Global Convergence
One SCC. Let us first consider a network that is strongly connected. In this case, the corresponding Jacobian matrix is irreducible. As all the off-diagonal elements are nonnegative, this implies that the dynamics of Eq. (2) is strongly monotone [18,19]. This allows us to apply the Hirsch convergence theorem, which gives that the system converges to a stable stationary solution if started from any initial condition (except for a set with zero measure) [18]. Using the Lyapunov function defined in the previous sections, we can further state that there is exactly one stable stationary solution: each such solution x Ã is attractive at least for initial conditions x §x Ã , which means that multiple stable stationary solutions would have overlapping basins of attraction.
The Hirsch theorem is valid only for strongly connected systems; now we present an alternative proof of the existence and attractive properties of the stable stationary solution of the system in Eq. (2) which can be readily generalized for systems which consist of multiple SCCs.
For simplicity we will refer to the right-hand-side of Eq.
(2) as f (with components f i : _ x x i ). Notice that in the absence of an external signal (s i~0 ), the origin is a stationary solution. If it is linearly stable, then the Lyapunov function (8) proves that it is globally attractive. Thus in the further analysis we will consider either the case when the origin is unstable (which we will refer to as case (i)) or that when there is an external signal (case (ii), note that now f i w0 for at least one i component). In these cases, we will first show that there is a nonzero stable stationary solution x Ã which can be connected to the origin with a path where fw0. Next we can prove that for each x initial condition which is not itself an unstable stationary solution, the system will reach the x §x Ã region where the Lyapunov function presented earlier is valid.
The existence of a stable stationary solution. We denote the Jacobian matrix at the origin by J 0 . Multiplying the Jacobian matrix with a unit vector e gives that how the components of f change, if we move in the direction pointed by e (i.e. f(ee)~f(0)zeJ 0 ezO(e 2 )).
In case (i), all derivatives in the origin are zero and the Jacobi matrix there has a real positive eigenvalue l with a corresponding positive eigenvector u. This implies that in the direction of u, all components of f are increasing; as all f i are continuous and infinitely differentiable, it follows that there is a set in the neighborhood of the origin, where fw0.
In case (ii), at least one derivative f i is positive at the origin (the others may be either positive or zero). In this case, we only need to find a direction where the zero components are increasing. This means that we need to find a positive vector v for which the vector w:J 0 v has w i w0 for values of i, where f i (0)~0 (the other components can be negative or zero in this case). Considering that J 0 is a Metzler matrix, this is always possible.
Thus in both cases, we get that there is a set P in the neighborhood of the origin, where all derivatives are positive: We further require that P is connected and that the origin is an accumulation point of P (0 [ P P). Now let us consider the closure of P, P P, which contains points, where f §0. Note that both P and P P are bounded: f i (x)v0 if x i~1 . Let Q be the set of points in P P which are ''farthest away'' from the origin: Now we can easily prove that every point in Q is a stable stationary solution. For this, let us consider a point x[Q. From the definition of Q, it follows, that there will be some points v[P which are close to x. All these points must satisfy DvDƒDxD. Now let us assume that x is not a stable stationary solution. This means that x is either an unstable stationary solution or a point where some components of f are positive, and some (or possibly none) are zero. In both cases, we will reach a contradiction by proving that there are some points y[P for which DyDwDxD. Using that the Jacobi matrix in x is a Metzler matrix and employing similar arguments as for J 0 , we can find a positive vector u which points in a direction where again all components of f become positive. A more precise formulation of this is the following. Without the loss of generality, we can assume that the first m components of f(x) are zero, and the rest are positive (the order of the components is arbitrary). Of course, m can be the total number of components in the case when x is assumed to be an unstable stationary solution. Using the very same arguments which we employed for the construction of the original P set, we have that Auw0, such that J(x)u ð Þ i w0 for all iƒm. Using that f is a continuous and differentiable function of its variables, we have that Aaw0 such that f(xzbu)w0,Vb[(0,a). Now we have to employ the fact that x[Q. This means that there are points in P which are close to x: Vew0,Av such that f(v)w0, DvDvDxD and Dx{vDve (of course, v[P). Furthermore, again using that f is continuous and differentiable, we have, that for small values of e, v can arise as a linear approximation: Adw0 and Aw such that J(x)w ð Þ i w0,Viƒm, and also, Vc[(0,d) we have DxzcwDvDxD and f(xzcw)w0 (i.e. we consider a line pointing from x into P).
Using these u and w vectors, we can now prove that Ab 1 w0 such that xzb 1 u:y[P. As uw0, we will have ywx, resulting in a contradiction. We have already seen that f (y)w0 for the appropriate choice of b 1 ; as we require P to be connected, we need to show that y can be connected by a continuous curve to some other points known to be in P. A such curve can be defined as z(t)~xzb 1 z 1 (t), where.
Notice that z(0)[P, for small enough b 1 . The derivatives along this curve will be: Considering the second term, we have that J(x)z 1 (t) ð Þ i w0, Vt[½0,1, and Viƒm (these components are just the sum of two positive numbers). This implies that Ab 1 w0 such that f(z(t))w0, Vt[½0,1. As z(0)[P, we have that z(t)[P, Vt[½0,1 also. Applying this to y:z(1) and using that DyDwDxD, we have reached a contradiction: if x[Q, no such y could exist (this was the definition of Q). As our basic assumption was that x is not a stable stationary solution, it follows that Vx[Q must be a stable stationary solution. Since P P is bounded and nonempty, also Q is nonempty, which means that there is at least one stable stationary solution. Applying the Lyapunov function presented previously, we get that there is exactly one stable stationary solution.
Convergence. Considering the P P set defined above, we have seen that it contains two special points of interest: the origin and the stable stationary solution x Ã of the system. (both of these are accumulation points of the set P). Using that the derivatives in Eq.
(2) are all continuous and infinitely differentiable functions, it follows that we can construct a curve r(t), connecting the origin and x Ã which is contained in the set P: There are of course infinitely many such curves, but any one of them is suitable. We will refer to an arbitrary such curve as r(t).
For every point on that curve, we will consider the following Ndimensional rectangle (if x has N components): For any of these R(t) rectangles in the 0vtv1 range, we have that the f derivative vectors point toward the ''inside'', in which the stable stationary solution resides. If any coordinate is x i~1 , its derivative will be negative. On the other hand, if some coordinate is x i~ri (t), the corresponding derivative will be positive, since f(r(t))w0 and the Jacobi matrix at any point has nonnegative offdiagonal elements. This implies that the solution of Eq. (2) is a flow, which heads into the x §x Ã region, meaning that the system started from any initial condition will eventually reach this region, where the Lyapunov function defined earlier guarantees its convergence to x Ã .
Further remarks. The above proof about convergence holds for initial conditions where xw0. It is easy to see that our arguments are also valid for initial conditions where x §0 given that either there is an external signal or x=0 (we exclude the origin which is an unstable stationary solution if there is no external signal). In this case we have that at least one derivative f i is positive. This means that there is a time t 1 such that the corresponding coordinate x i will be positive. That implies that for all j, where a ji w0, the f j components will also be positive, meaning that at some time t 2 wt 1 the x j coordinates affected by x i will be positive. Using that the network is strongly connected, and the derivatives are continuous and infinitely differentiable, we get that after some time t, all components will be positive. Since we are in the xw0 region now, the previous proof can be applied in this case too.
We also note that we have excluded the possibility when the largest eigenvalue of the Jacobian matrix at the origin is exactly zero. In this case our previous arguments do not hold. We think that in a real biological system this is not a concern; Eqs. (2) or (3) are valid only if the time scale of their solutions matches the time scale arising in other biological processes in the cell (i.e. a model based on ordinary differential equations is clearly unable to model the system if l (0) t%1, where l (0) is the eigenvalue of the Jacobian matrix at a stationary solution and t is the relevant timescale in the system).
More SCCs. We can inductively extend the proof presented in the previous section to a network which is not strongly connected, but is made up of an arbitrary combination of SCCs. Let us consider a network which is made up of two subnetworks, C 1 and C 2 , with links only in the C 1 ?C 2 direction between them.
We need to prove that if our results hold for both C 1 and C 2 then they will also hold for the whole network. We can order the variables in a way that C 1 and C 2 form blocks in x. We will denote the variables belonging to each component with subscripts (x 1 and x 2 for the independent variables and f 1 and f 2 for the derivatives). Let us assume that our previous results hold for both components: for i~1,2, the system C i has a stable stationary solution in the presence of any stationary inputs. Also if this stationary solution is not the origin, then it can be connected to the origin with a path, where all derivatives are positive.
Let x Ã 1 be the stable stationary solution of system C 1 . Using x Ã 1 as constant inputs on system C 2 , we can determine the x Ã 2 stable stationary solution of C 2 . Now, the vector x Ã~( x Ã 1 ,x Ã 2 ) will be the stable stationary solution for the whole system.
Let us consider the Jacobian matrix in the origin, J 0 , and also the submatrices corresponding to the two subsystem, J (1) 0 and J (2) 0 . Let l (i) denote the largest eigenvalue of J (i) 0 . We have three separate cases: (1) Case (1) is covered by the Lyapunov function, so we have to give the proof only for cases (2)- (3).
Case (2) The component C 1 relaxes to its stationary solution. If we set x 1~x Ã 1 :0, then the proof presented in the previous section is valid for component C 2 . If we now set x 1 w0, then all components of f 2 increase. Thus, the arguments presented for convergence in the previous section remain valid; the C 2 component will reach the x 2 §x Ã 2 region. After that, for the whole system, we will have x §x Ã , which again allows us to employ our Lyapunov-function, getting that the whole system converges to x Ã .
Case (3) The proof of the global attractive properties is similar to the case of a strongly connected network, we only need to change some properties of the r(t) curve connecting the origin and the stable stationary solution. First we define the set P 1 for the system C 1 in a similar way as the set P in the previous section (see Eq. (17)). We also define the set P 2 for system C 2 in this way with the restriction that we consider the system C 2 with constant inputs x Ã 1 . We consider the curve r(t)~(r 1 (t),r 2 (t)), and require that: Again, we can use the continuously shrinking rectangles R(t) to prove that the system will reach the x §x Ã region and then converge to x Ã .
This way we get that if our previous results hold for systems C 1 and C 2 , then they also hold for the combined system C 1 ?C 2 . A network with more than two SCCs can always be built with the addition of one SCC at a time. Thus, for any network, we can start from an individual SCCs, and iteratively add each SCC, applying the arguments presented in this section at each step, getting that our results apply for the whole network. Considering the simple example network presented in Fig. 3., this means that we can prove our results for the whole network starting with SCC1 and consecutively adding SCC2, SCC3, SCC4 and SCC5.

Discussion
We have shown, that the model suggested for the general treatment of signaling networks in [10] (Eq. (2)) always has one stationary solution, which is globally attractive (any initial condition relaxes to it), which means that the only relevant qualitative property of that model is that whether the off-state or an autoactivated state is the stable stationary solution in the absence of an external signal. This gives an example for a class of nonlinear systems which show only very simple behavior, and also an example for nonlinear systems where a Lyapunov function can be employed proving the attractive properties of a fixed point. This behavior is a consequence of the parameter space being bounded and the monotonic nature of Eq. (2), which gives a very simple flow as a solution. We have also shown that the straightforward generalization of the system to include inhibitory interactions (Eq. (3)) does not change the behavior in the case when the off-state is stable, thus the stability of the off-state can be determined by analyzing the network with only the positive interactions present.
While a more general model can have complex solutions, the observation that the stability of the off-state is determined by the behavior of the dynamically positive cycles in the network applies to a wider range of possible models. Counting such cycles in the Signalink database [4] of the human signaling network and also in randomly generated networks, we found that the real signaling network has significantly few dynamically positive cycles, which agrees with our expectations. This result supplements the findings of Kartal and Ebenhöh [12], who gained similar results for network data including only positive interactions.
Based on our findings, we expect a system transmitting an external signal to have a stable off-state which it returns to if the external system has been shut off. This should be the correct behavior in many cases. If, for some reason, this behavior changes the cell will ''think'' that there is an external signal when only the signal transmitting network is stuck in a faulty autoactivated state [12]. This could have drastic effects on the cell and cause it to cease carrying out its original purposes. The stability properties of the offstate could change following a change in the parameters of the system, caused by either the mutation of the proteins involved or the change in the protein concentrations. In the simple model considered here (Eq. (2) and (3)), this means a change in the matrix elements a ij f or c ij . Computing the Jacobian matrix at the off-state and using the Perron-Frobenius theorem and the Collatz-Wielandt formula [17] we can see that the increase in the activation rate constants or the total protein concentrations causes an increase in the largest eigenvalue, and an increase in the deactivation rate constants causes the largest eigenvalue to decrease, while the inhibitory interactions c ij do not affect the linear stability. Mutations in certain proteins can thus lead to either the off-state of the network losing its stability giving a cell a constant stimulation as if it was constantly under an external signal, or to the dampening of the signal, lessening its effects on the cell (or even to the silencing of a pathway if a link is completely taken out). In our simple model a change in one of the parameter values can be compensated by changing some of the other parameters (e.g. the concentration of the phosphatases present). In a real network, more sophisticated methods will be needed, but we believe that our results concerning the importance of dynamically positive cycles can lead to a better understanding of intracellular signaling networks.