Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS Comput Biol**|**v.4(3); 2008 March**|**PMC2265423

Formats

Article sections

Authors

Related links

PLoS Comput Biol. 2008 March; 4(3): e1000041.

Published online 2008 March 21. doi: 10.1371/journal.pcbi.1000041

PMCID: PMC2265423

Herbert M. Sauro, Academic Editor^{}

University of Washington, United States of America

* E-mail: rf.srnc.nlni@erhclupeS.erdnaxelA-seuqcaJ

Conceived and designed the experiments: ACV J-AS. Analyzed the data: ACV J-AS. Wrote the paper: ACV J-AS SDM.

Received 2007 September 21; Accepted 2008 February 20.

Copyright Ventura et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

This article has been cited by other articles in PMC.

Cycles involving covalent modification of proteins are key components of the intracellular signaling machinery. Each cycle is comprised of two interconvertable forms of a particular protein. A classic signaling pathway is structured by a chain or cascade of basic cycle units in such a way that the activated protein in one cycle promotes the activation of the next protein in the chain, and so on. Starting from a mechanistic kinetic description and using a careful perturbation analysis, we have derived, to our knowledge for the first time, a consistent approximation of the chain with one variable per cycle. The model we derive is distinct from the one that has been in use in the literature for several years, which is a phenomenological extension of the Goldbeter-Koshland biochemical switch. Even though much has been done regarding the mathematical modeling of these systems, our contribution fills a gap between existing models and, in doing so, we have unveiled critical new properties of this type of signaling cascades. A key feature of our new model is that a negative feedback emerges naturally, exerted between each cycle and its predecessor. Due to this negative feedback, the system displays damped temporal oscillations under constant stimulation and, most important, propagates perturbations both forwards and backwards. This last attribute challenges the widespread notion of unidirectionality in signaling cascades. Concrete examples of applications to MAPK cascades are discussed. All these properties are shared by the complete mechanistic description and our simplified model, but not by previously derived phenomenological models of signaling cascades.

Cellular signaling is carried out by a complex network of interactions. A structure that is found commonly in signaling pathways is a sequence of on-off cycles between two states of the same protein, referred to as a cascade. By analyzing and reducing the basic kinetic equations of this system, we have constructed a new mathematical model of an intracellular signaling cascade. It is widely accepted that information travels both outside-in and inside-out in signaling pathways. Conversely, cascades, even while being main components of those pathways, have been so far understood as structures where signal transmission occurs in a manner analogous to a domino effect: the information flows in only one direction. Adding explicit connections linking a particular level with an upstream location has been the way bidirectional propagation has been explained so far. In other words, up to now, unidirectional cascades would allow bidirectional propagation only when embedded in more complicated circuits. The proposed model shows that a cascade can naturally exhibit bidirectional propagation without invoking extra re-wiring. This result inspires novel interpretations of experimental data; since signaling pathways are usually reconstructed from such data, this outcome could have far-reaching implications in the understanding of cell signaling.

Covalent modification cycles are one of the major intracellular signaling mechanisms, both in prokaryotic and eukaryotic organisms [1]. Complex signaling occurs through networks of signaling pathways made up of chains or cascades of such cycles, in which the activated protein in one cycle promotes the activation of the protein in the next link of the chain. In this way, an input signal injected at one end of the pathway can propagate traveling through its building-blocks to elicit one or more effects at a downstream location.

Examples of covalent modification are methylation-demethylation, activation-inactivation of GTP-binding proteins and, probably the most studied process, phosphorylation-dephosphorylation (PD) [1],[2]. In such cycles, a signaling protein is activated by the addition of a chemical group and inactivated by its removal. This protein is modified in turn by two opposing enzymes, such as a kinase and a phosphatase for PD cycles. In the absence of external stimulation, a cycle exists in a steady state in which the activation and inactivation reactions are balanced. External stimuli that produces a change in the activity of the converting enzymes, shifts the activation state of the target protein, creating a departure from steady state which can propagate through the cascade.

The advantages of these cascades in signal transduction are multiple and the conservation of their basic structure throughout evolution, suggests their usefulness. A reaction cascade may amplify a weak signal, it may accelerate the speed of signaling, can steepen the profile of a graded input as it is propagated, filter out noise in signal reception, introduce time delay, and allow alternative entry points for differential regulation [3]–[5].

Intracellular signaling through cascades of biochemical reactions has been the subject of a great number of studies (e.g., [2],[6] for reviews). Theoretical investigations have been motivated by the increased need for developing an abstract framework to understand the vast amounts of experimental data now available. This whole field of research is further motivated by the hope of characterizing pathways that are deregulated in diseases such as cancer and to define targets to combat these diseases [7].

Since the stimuli a cell receives are varied and complex, cascades do not operate in isolation, but rather the integration of stimuli depends on crosstalk between pathways. Another crucial property of signaling cascades is their ability to integrate information by transmitting the effects downstream and also feedback upstream. In spite of a few decades of intense work on signaling cascades, no models have ever been built that exhibit crosstalk with backwards and forwards transmission of a lateral input from another cascade, except when *ad hoc* feedback is explicitly added to the cascade model. Our model, built from first principles, naturally exhibits these characteristics and therefore inspires novel interpretations of experimental data.

A well studied example of a cascade of activation-inactivation cycles is the cascade of protein kinases. In this case, the basic signaling unit is a PD cycle, whose activating kinase is the phosphorylated protein of the previous cycle. Many proteins contain several phosphorylation sites, allowing for great versatility of regulation. Such is the case, for example, for the mitogen-activated protein kinase (MAPK) cascade, which is widely involved in eukaryotic signal transduction [3], [8]–[10]. For the sake of simplicity, in this article we will mostly consider cascades composed of simple, 2-state activation-inactivation cycles. However, the equations corresponding to the MAPK cascade are also derived and some of their properties compared with those of the simpler cascades. Even though our results are valid in general, for covalent modification cycles, we will employ the nomenclature associated with PD cycles, *i.e*. the converting enzymes will be referred to as kinase/phosphatase.

The focus of our study is to refine the mathematical modeling of cascades of covalent modification cycles, such us the one depicted in Figure 1. Several mathematical descriptions have been developed to describe such cascades using ordinary differential equations. Typically, those descriptions are built up starting with a model for a single cycle, which is then phenomenologically incorporated into a cascade of cycles. A well known model for describing the single cycle was introduced by the pioneering work of Goldbeter and Koshland (GK) [11]. The GK model considers the concentration of the target protein to be in large excess over those of the converting enzymes, thereby reducing the description to a single equation per cycle. The model obtained in this way was then phenomenologically extended to a cascade of individual GK cycles. Here, by the designation “phenomenological” we mean that, in the cascade, the forward coupling between the GK cycles is chosen as simply as possible, but not strictly deduced from first principles. This phenomenological framework extension of the GK model will be denoted as the *GK-like* model. The GK-like model has been used by several authors to describe the dynamics of signal transduction [9], [12]–[16]. For particular limiting cases, the GK-like model can be simplified further, which results in a model where the inter-converting reactions follow linear rate laws with first-order rate constants. This description was studied in several key papers [17]–[19], and we will refer to it as the *linear rates* model.

The concept of a “cascade” in the study of transduction pathways is appealing because of its modular structure. What is especially appealing is the possibility of defining the cascade state by only one variable per module. As mentioned above, since the building blocks of the GK-like model are the well-studied GK cycles, they involve only one equation per cycle. A different approach however, is to deal with the dynamics of the cascade of Figure 1 by considering the complete set of biochemical reactions and by writing the corresponding equations without any upfront approximations. This was accomplished, for example, for the case of the MAPK cascade [8]. We will refer to this approach as the *mechanistic* model. For the purposes of this paper, we will consider that the mechanistic model represents a complete description of the system under study (event though we recognize that, in reality, it is not a hypothesis-free model).

In this article, starting from the mechanistic description of a cascade composed of an arbitrary number of cycles, we derive a consistent approximation under which the cascade is described with one variable per cycle. It turns out that in this derivation, referred to as a *reduced mechanistic* description, the phenomenological GK-like model is not recovered. At first sight, our new approximation differs slightly from the previously derived description for signaling cascades. However, it involves qualitatively different dynamics from the GK-like model, yet it is in very good agreement with the complete mechanistic description when the approximation conditions are fulfilled.

The main difference between our simplified mechanistic description and the phenomenological one is the appearance of an intrinsic feedback from each unit to the preceding one, caused by the fact that in each cycle there is sequestration of part of the activated protein of the previous step. The new description of the cascade predicts the existence of damped oscillations along the chain, a phenomenon that cannot be observed using the previous phenomenological description. Interestingly, a corollary of our model is that if a particular unit in the middle of the chain receives an input–a common event, given the high degree of crosstalk between signaling pathways–then our reduced mechanistic description predicts that this perturbation is able to travel both forwards and backwards. This “bicistronic” propagation, which may be critical for effective eukaryotic signaling, is not possible within the GK-like description either. Our model provides a suitable framework for future experiments that investigate crosstalk and bicistronic propagation of signals.

We consider a cascade of biochemical cycles, as illustrated on Figure 1, in which the two variables *Y _{i}* and

Since processes involved in the control of production of new proteins proceed at a much slower timescale than processes that chemically modify existing proteins, the total quantity *Y _{iT}* is considered to be constant in time and the variables [

In an ideal situation, the inter-conversion of the *i ^{th}* protein can be described by the following reactions:

(1)

where *C _{i}* and

Working in the framework of the mechanistic model offers the advantage that no mathematical approximations are needed (even though, overall, this is obviously not a hypothesis-free model), and this could be the optimal choice for comparing experimental data with numerical simulations of the model. This option was taken, for example, in the context of the MAPK cascade [8].

On the other hand, more complicated models, although in principle more realistic, are also less amenable to developing insights into the transduction pathways. It is appealing then, to find out under which set of hypotheses can the mechanistic model be approximated by a simpler one, for example a model with a single variable per cycle. Arriving to such reduced description is the main contribution of the present paper. Before providing that description, we briefly review some approaches followed in the literature to study signaling cascades by means of only one equation per cycle (see Text S2 for a summary).

One possible simplification of the chemical reactions in Equation 1 is to neglect the formation of the complexes *C _{i}* and

(2)

with the definitions *y _{i}*=[

The linear rates model does not account for the fact that the transformations from *Y _{i}* into

(3)

where *y _{i}*+

(4)

We will refer to this generalization of the model of Goldbeter and Koshland as the GK-like model. Let us note that in the case where the coefficients *K _{i}* and

Equation 3 was first derived by Goldbeter and Koshland, to study the so-called property of zero-order ultrasensitivity. This means that when the *K*′s are small (e.g., of the order of 10^{−}^{2}) the cycle behaves like a switch where the steady state for *y ^{*}* passes abruptly from its lowest to the highest value as a function of the ratio

The cascade extensions of the Goldbeter-Koshland model have been extensively used in several important articles [9], [12]–[16] covering different contexts, often adding a single negative feedback loop extending from the last unit of the chain to the first one. It has been argued however, that the hypotheses leading to the chain equations (Equation 4) are questionable [3],[20]. Blüthgen et al [3] claim that, in several cases, current experimental data do not support neglecting the enzyme-substrate complexes from the conservation equation. If the affinity of kinase *Y _{i}^{*}* for the protein

The importance of sequestration-based feedback in signaling cascades is thoroughly analyzed in the recent work by Legewie et al [21], where a positive feedback mechanism that emerges from sequestration effects is shown to bring about bistability in the cascade. In that study, sequestration is caused by stable heterodimers formed by the non-phosphorylated protein *Y _{i}* and the next substrate

In Text S1 we derive in detail the new class of model equations obtained as an approximation of the mechanistic model. The goal of our approach is to reduce the number of variables in the complete system by bringing into play hypothesis that allow us to use the quasi-steady state approximation. Three key dimensionless parameters are defined to facilitate the analysis:

(5)

*ε _{i}* and

Using a standard singular perturbation analysis, we have found that the state of each biochemical cycle can be described by a single variable defined as *x _{i}=y_{i}^{*}+c_{i}*

(6)

with the following conservation equation from which *y*_{i} has to be extracted:

(7)

*x*_{0}*=S* is the normalized input signal and *y _{n}*

The reduced system given by Equations 6–7 seems to be, in principle, equivalent to the GK-like model given by Equation 4. However, two main features make it significantly different. First, in our novel system, termed the reduced mechanistic model, the conservation equation depends on the variable of the previous cycle. Second and more interesting, the denominator of the negative term in Equation 6 is now a function of the next variable *y _{i}*

In the limit *η _{i}~ε_{i}«*1, one retrieves the simple conservation law

In addition to the limit *η _{i}~ε_{i}«*1, our perturbation scheme encompasses situations where the total protein does not necessarily increase along the cascade. We then allow

Finally, we notice that our description enables a reduction of the cascade equations with mixed hypothesis concerning the enzymatic reactions. For example, we could have *µ*_{1}*~ε*_{1} and *η*_{1}*~*1 for the first cycle, *µ*_{2}*~*1 and *η _{2}~ε_{2}* for the second cycle, etc. Or even

In Text S3, we present the extension of the reduced mechanistic model for a cascade involving double-phosphorylation. Notwithstanding that these equations are more complicated than Equations 6–7, the distinguishing feature is maintained: each level in the cascade is subject to influence from the following level which, in the appropriate *x _{i}* variable, can be identified as a negative feedback. In the current study we analyze mostly static properties of these more complicated equations and compare them to those of Equations 6–7, while a more exhaustive characterization will be presented in a future article.

In this section we report on dynamic and static properties of the new chain equations (Equations 6–7), when studied by numerical simulations, and compare them to those of previous cascade models. We will consider both short (*n*=3) and long chains (*n*=10, or 15), respectively. In all the figures we plot *y _{i}^{*}*, the level of active protein, obtained from

The homogeneity assumption implies that *V _{i}*

In Figure 2 we present an initial exploration of the dynamics that the reduced mechanistic model is able to display and how well it approximates the mechanistic model. As an example, a 10-unit chain was considered. The temporal evolution of the variable describing the first unit, *y*_{1}* ^{*}*, is plotted. For each choice of

We have also computed the errors for less extreme conditions, such as 1) *ε*=0.1, *η*=*µ*=0.5; 2) *ε*=0.1, *η*=0.5, and *µ*=1; and 3) *ε*=0.1, *η*=1, and *µ*=0.5 (notice that *ηµ*~2.5 *ε* for 1) and *ηµ*~5 *ε* for both 2) and 3)). The respective errors are 4.3%, 3.5%, and 5% (data not shown). The errors in the prediction of the steady states are lower than 0.0001% in all the mentioned cases, indicating the high accuracy of the reduced mechanistic model to study static features of the cascade. This property is due to the conservation equation, Equation 7, taking into account the first correction in *ε _{i}* (see Text S1).

Interestingly, the temporal evolution of activated protein depicted in Figure 2 exhibits damped oscillations. This feature is displayed by both the mechanistic model and the novel reduced description introduced here. However, such behavior is not attainable within the other available descriptions for signaling cascades models, *i.e*. Equations 2 and 4. Our simplified new description reveals this attribute of the complete (mechanistic) model, that has remained (to our knowledge) hidden until now.

The existence of damped oscillations has been corroborated by a numerical study of stability of the steady state of Equation 6, which is a stable focus. The spectrum of the Jacobian matrix of this system computed at steady-state indeed possesses several eigenvalues with nonzero imaginary parts. The real parts however, are always negative (as observed, e.g., by continuation in *S* parameter) and therefore we cannot obtain sustained oscillations in this chain. A more detailed mathematical study of the spectrum of stability of the chain is beyond the scope of this paper and will be the object of future work. Damped oscillations are not possible without a negative feedback between the cycles [24] and thus reflects the new feature of our model of cascades.

Figure 3 contains a representative characterization of the model's temporal dynamics. To simplify the description, we consider two control parameters: *V*=*µη/ε* and *η*, while the other parameters are set as *ε*=0.01, *K*=0.01 and *K′*=*K/µ*. The input stimulation is turned on at time 0 from *S*=0 to *S*=1. Figure 3A displays the parameter space *η*−*V*. In every panel of Figure 3B, the temporal evolution *y _{i}^{*}* for three of the units in the chain are plotted. In some of the plots, the

Another interesting and surprising feature revealed by our new model is that, even though the parameters are homogeneous through the chain, the steady states of variables *y _{i}^{*}* do not always exhibit a monotonous trend with respect to the unit index

Let us consider a chain that is in equilibrium, *i.e.*, a cascade in which every unit has achieved steady-state. We then perturb a single variable *x _{i}*, as indicated in Figure 4. The nature of our new description, that couples each unit with both the previous and the following one, makes it possible to transduce the localized perturbation in both directions, forward and backward. Figure 4A and 4B, which correspond to parameters A

We now study the time-independent features of our model (Equations 6–7), and compare them with those of the system described by Equation 4 for the same set of parameters. The stimulus-response curve is defined, as customary, by the steady-state values of the variables as a function of the input stimulus *S* (recall that in Equation 6, *x*_{0}=*S*, which is assumed to be constant in time). Figure 5 shows the stimulus-response curves obtained with Equation 6 (filled lines) for condition A_{1} from Figure 3A, except for the dotted line that was calculated with the GK-like model (Equation 4). We observe that with the parameters so assigned, the computations performed with the reduced model deviate by less than 0.0001% from the (non approximated) mechanistic model. The results were obtained with a chain of three units, as an example.

In Figure 5, the variables *y _{i}^{*}* display sigmoidal responses. The low level of activation for units 1 and 2 is due to the fact that the proteins are indeed partially sequestered in the enzyme-substrate complexes. However,

In this section we apply the reduced mechanistic model to a well-known signaling pathway, the mitogen-activated protein kinase (MAPK) one [3], [8]–[10],[25],[26]. We first base our description on a particular published set of parameters for this pathway [8]. Importantly, the results obtained are not qualitatively modified by variations of the selected values in the ranges suggested in the literature [8]. Moreover, they are not modified by choosing different sets of parameters [9],[25],[26], as described in the Text S6.

It is well know that the MAPK cascade consists of three levels, the second and the third ones involving a double-phosphorylation mechanism. In this section we consider both the MAPK cascade and a simpler case, a 3-unit chain where each unit is a 2-state cycle.

Starting with the published set of parameters (see [8] and also [10], for a summary), we have computed the parameters involved in the reduced mechanistic description and listed them in Table 1. As described in Text S6, there are some extra parameters for the case involving double-phosphorylation, that are designated *ν*, *K ^{*}*, and

According to Table 1, the conditions under which the reduced model is valid are only partially satisfied, *ηµ*~*ε* for the first unit but *ηµ*~10 *ε* for the second and third ones. Even for these conditions and since the focus of this section is in steady states, the reduced mechanistic model provides a description that is in excellent agreement with the complete mechanistic one.

In Figure 6 we plot the normalized stimulus-response curves for a 3-unit chain, either with single-phosphorylation in all the units (A) or with single-phosphorylation in unit 1 and double-phosphorylation in units 2 and 3 (B), *i.e.*, the case corresponding to the MAPK cascade. Both cases are characterized by the parameters in Table 1. The input stimulus was taken to be the concentration of *E*_{1T} , the total amount of kinase for the first unit in the cascade (corresponding to MAPK kinase kinase in B). *E*_{1T}, related to the parameter *S* we have used as input in the previous section, was varied over a wide range. The outcomes were obtained by both the complete mechanistic and the reduced mechanistic models and the results are indistinguishable for the scales of the figure (black, blue, and red filled lines for *y*_{1}^{*}, y_{2}* ^{*}*, and y

In order to compare the steepness in the responses, we have computed the apparent Hill coefficient *n _{H}* ([8]) for each curve, as indicated in the legend. As expected,

We have mentioned that the good agreement between the mechanistic and reduced mechanistic descriptions regarding the prediction of steady states is due to the conservation law, Equation 7, taking into account the first correction in *ε _{i}* (see Text S1). If that correction is not considered, differences could appear in the steady states predicted by the mechanistic model and the reduced mechanistic one. However, and for the parameters in Figure 6, the predicted values of

In Figure 6B, the mechanistic and reduced mechanistic models' outcomes and corresponding Hill coefficients recover published results [8]. Comparing figures A and B, we also confirm that the chain involving double-phosphorylation responds in a steeper manner than the one with only single-phosphorylation, as expected from previous work [27].

In Figure 7 we show the outcome of stimulating the 3-unit chain as indicated in the schemes close to each panel: the input stimulus to the cascade was taken to be the concentration of *E′*_{3T}, the total amount of phosphatase for the last unit in the cascade (corresponding to MAPK in B). *E′*_{3T} was varied over its suggested range of variation [8]. Increasing the amount of phosphatase produces a decrease in the response curve *y*_{3}* ^{*}* (red filled line), as expected. Interestingly, our new reduced model (Equations 6–7), as well as the complete mechanistic description, predict that this perturbation on the third level of the chain is propagated backwards: the variation in

The insets in both figures indicate that is not necessary to vary parameter *E′*_{3T} over a wide range to observe this property, rather it is clearly seen by changing it only by a factor 5 around its suggested concentration (0.12 *µM*), where a 20% variation in *y*_{2}* ^{*}* is observed, a value that is high enough to be detected experimentally (meaning that it is most likely not contained within the error of the experiment). Due to the parameters characterizing this particular pathway, the effect is not propagated to

In Text S6 we extend the results in this section concerning “reverse” stimulus-response curves, for different sets of published parameters on the MAPK cascade.

A modular response analysis (MRA) [28] was applied to determine the network architecture of the cascade in the context of the new model equations (Equations 6-7). MRA has recently been proposed as a tool to characterize the interactions between “modules” in a cellular regulatory network, having the advantage of allowing direct experimental implementation.

As a matter of fact, the negative sign of the Jacobian element *x* ˙_{i}/x_{i}_{+1} indicates that the (*i*+1)* ^{th}* level of the cascade exerts a negative effect in what concerns variable

As a result of applying MRA, a matrix of local response coefficients **r** is obtained. An element *r _{ij}* in this matrix describes how the state of the variable associated with module

Indeed, if the rate of change of variable *x _{i}* is denoted by the function

(8)

meaning that *r _{ij}* corresponds to a scaled version of the Jacobian matrix

Using notations and concepts from the literature [28], we apply the MRA method to a 3-unit cascade involving only single-phosphorylation and characterized by the parameters in Table 1 [8]. There are three modules in this network as described by Equations 6-7, each of them corresponding to the three successive levels in the cascade and characterized by a single variable *x _{i}*. Figure 8A contains the matrix of local response coefficients

The structure of matrix **r** is tridiagonal, meaning that the first level in the cascade does not directly influence the third one (*r*_{31}=0), and viceversa (*r*_{13}=0). Coefficients *r*_{21} and *r*_{32} are positive, representing the positive effect of each level in the cascade to the subsequent one. Interestingly, *r*_{12} and *r*_{23} are both negative, indicating an inhibitory effect from unit (*i*+1) to unit *i*. The resulting connections between the units in the cascade are summarized in the scheme in Figure 8A.

To understand these results in more depth, we have studied how the coefficients in the matrix in Figure 8A depend on the parameters characterizing the cascade. For example, Figure 8B shows coefficients *r*_{21}, *r*_{32}, *r*_{12}, and *r*_{23} versus the parameter *E*_{1T}. *r*_{31} and *r*_{13} are zero (data not shown), *r*_{21} and *r*_{32} are positive, and *r*_{12} and *r*_{23} are negative, throughout the range where *E*_{1T} was varied. Depending on the value of *E*_{1T}, each of the nonzero *r _{ij}* could be less or greater than 1 and the relative strength of the backward and forward couplings for a given pair of modules, e.g. |

Studies like the one in Figure 8B help us understand and also predict the degree of backwards coupling as a function of the parameters in the model. One utility of this work is as a starting point of a more systematic study on how to enhance or attenuate that coupling in the cascade, the subject of our ongoing work.

Interestingly, the interaction map characterizing the connectivities between variables *x _{i}* (matrix

Experimental data concerning the application of MRA to the MAPK cascade are now available in the literature [30], showing non zero *r*_{21} and *r*_{32} coefficients (and also non zero *r*_{31} and *r*_{13} coefficients). The interpretation of non zero *r*_{31} and *r*_{13} was proposed in terms of the usual “explicit” positive or negative feedbacks which are sometimes considered in models of signaling pathways [9],[12],[13]. From this perspective, the explanation for the non zero *r*_{12} and *r*_{23} coefficients was, at least regarding *r*_{23} and based on experimental evidence, that not only is *y*_{2}* ^{*}* able to phosphorylate

The main contribution of this work is to propose a new one-variable per cycle model for signaling cascades of covalent modification cycles, consistent with a mechanistic complete description. Our model reveals new and biologically relevant properties of such cascades. These properties are characterized completely for the case of single-phosphorylated cascades. Furthermore, single and doubly-phosphorylated cases are compared regarding their stimulus-response curves, while a more exhaustive characterization of the scheme involving double phosphorylation will be presented in a future article.

The scheme in Figure 1, which has been employed by many groups, is suggestive of the concept of a “cascade”. From a systemic point of view, a cascade is a system composed of units, the output of which is successively an input to the next unit. Based on this structure, powerful concepts from control theory can be applied successfully to the study of signaling cascades [14]. Although these concepts have proven its utility in many contexts, this kind of schematic representation implicitly conveys the idea that a signaling cascade is only a feed-forward chain in which signal transmission is analogous to a domino effect [32],[33]. Our study sheds a different light on this system, showing that this schematic representation can be misleading, since it turns out that each unit is actually coupled not only to the following one but also to the previous one, and interesting dynamics can arise from these interactions.

Our initial motivation for developing a new one-variable description of signaling cascades, was the following observation. The main assumption underlying the GK description of a single cycle is that the concentration of the target protein is in large excess compared to those of the converting enzymes. Holding the same assumption over a cascade of units would mean that the target proteins are in higher and higher concentration as the cascade progresses, since they act as the transforming enzyme for the following cycle. To our knowledge, this important issue has not been remarked upon in the literature, except for a brief comment in the work of Millat et. al. ([20], page 11).

In order to get more insight into this point, we have sought special limiting cases for which the mechanistic and the GK-like model are in good agreement. However, it turns out that the dynamics of the signaling cascade described by the mechanistic and the GK-like models cannot be compared consistently. The fundamental reason for this mismatch is that a careful perturbation analysis applied to the mechanistic model provides a different set of equations.

We note that in search for an adequate set of hypothesis leading from the mechanistic equations to the model given by Equations 4, we have studied an alternative scheme in which the modified protein *Y _{i}*

Our mathematical method relies on the standard quasi-steady state assumption (QSSA), which can be applied under well defined conditions to elicit a clear separation between the slow and fast dynamics of the mechanistic model. Under this standard QSSA framework, our analysis shows that a good slow variable for which evolution equations can be written is the sum of the free activated enzyme which is available in the *i ^{th}* cycle plus the amount of this protein which is captured by the next inter-converting cycle. The idea of working with a mixed variable

Even in the less extended QSSA framework, the conditions under which the model is valid are made clear. Under such conditions, our new model is indeed in perfect agreement with the complete mechanistic model (Figure 2). Those conditions are expressed in terms of three key parameters (Equation 5) we have defined to simplify the study. Even though the phenomenological equations, Equation 4, are appealing because of their simpler form and modular nature, we could not find any set of assumptions that would enable us to recover those descriptions. Our simplified model reveals properties of signaling cascades that were either hidden by the complex structure of the complete mechanistic model or lost in the simplified phenomenological descriptions.

It was stated that the reduced mechanistic model is valid whenever these two conditions are satisfied: *ε _{i}«*1 and

All the novel properties of a signaling cascade reported in this paper are linked, as previously mentioned, to the negative feedback from each unit to the previous one. This backward negative feedback can produce damped temporal oscillations in the chain, or it can create amplified “pathway” oscillations in the steady states of the cascade. Interestingly, it can also transduce a signal both forward and backwards. Given the multi-branched complex nature of many signal transduction pathways, this finding could have wide implications and can help focus further experimental investigation.

It has recently been reported that the 3-level MAPK cascade has autonomous oscillations without any kind of added explicit feedback [36]. Following a systematic numerical exploration of the corresponding mechanistic model [8], the authors provide a qualitative description of the mechanism responsible for these sustained oscillations. Their explanation strongly suggests the necessity of a bistable behavior at the second or third levels of the cascade, thus requiring double-phosphorylation at these stages [37]. Consistent with their findings, we have observed only damped oscillations in the dynamics of the single-phosphorylated cascade (Equations 6–7), which has been the main focus of the present work. Interestingly, preliminary numerical simulations of our reduced doubly-phosphorylated cascade model (Text S3), indicates that these autonomous oscillations are recovered in the simplified description.

The stimulus-response curves of the new model were also investigated (Figure 5). They have the usual sigmoidal shape characteristic of ultrasensitive responses; however, they exhibit lower steepness when compared with the output of the GK-like model. This result corroborates the conclusions stated in the work of Blüthgen et al. [3], where an analysis of the effect of sequestration was conducted. This effect is partially mitigated by double-phosphorylation (Figure 6), as expected from the literature [27].

To further characterize the new model within realistic conditions, we have studied it subject to different sets of published parameters corresponding to a well-known signaling pathway, such as the MAPK one (Figures 6 and and7,7, and Text S6). We have found that the ability of the model to transduce a signal both forward and backwards is widespread and that the effect is of enough magnitude to allow experimental verification.

Finally, we have applied a modular response analysis to determine the network architecture of the cascade described by the new model equations (Figure 8). This well-known approach enables not only to test the bidirectional structure of the cascade, but also to estimate the relative strength of the backward interaction.

In summary, our findings do not at all weaken the importance of previous models like the GK-like models and those with linear rates. To the contrary, the results of our model provide a different approach to deal with a simple one-variable per cycle model to describe signaling cascades. We hope that our contribution will help in the understanding of existing models for signaling cascades, will improve the description of available data, and will inspire both theoretical and experimental investigation.

All the ODEs were integrated in MATLAB 7 (Mathworks, Natick, MA). The stimulus-response curves were obtained using MATCONT, a MATLAB package for numerical bifurcation analysis of ODEs. The symbolic calculations were done using the Symbolic Math Toolbox in MATLAB.

Perturbation analysis of the mechanistic model.

(0.06 MB PDF)

Click here for additional data file.^{(57K, pdf)}

Available models for signaling cascades using one variable per unit.

(0.06 MB PDF)

Click here for additional data file.^{(61K, pdf)}

Cascades involving double phosphorylation.

(0.23 MB PDF)

Click here for additional data file.^{(224K, pdf)}

Comparison variables *x _{i}* and

(0.49 MB PDF)

Click here for additional data file.^{(481K, pdf)}

About the amplified “pathway” oscillations.

(0.05 MB PDF)

Click here for additional data file.^{(47K, pdf)}

More on “reverse” stimulus-response curves.

(0.53 MB PDF)

Click here for additional data file.^{(521K, pdf)}

A variation of the model with an intermediate step in the kinase activation.

(0.21 MB PDF)

Click here for additional data file.^{(203K, pdf)}

We would like to thank the editors and reviewers who helped us improve the article.

The authors have declared that no competing interests exist.

This work was supported by grants from the Department of Defense Breast Cancer Research Program (ACV), The Burroughs Wellcome Fund (SDM), the Breast Cancer Research Foundation (SDM) and the NIH (RO1 CA77612).

1. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, et al. Molecular biology of the cell. Fourth Edition. Garland Science 2001

2. Kholodenko BN. Cell-signaling dynamics in time and space. Nat Rev Mol Cell Biol. 2006;7:165–176. [PMC free article] [PubMed]

3. Blüthgen N, Bruggeman FJ, Legewie S, Herzel H, Westerhoff HV, et al. Effects of sequestration on signal transduction cascades. FEBS J. 2006;237:895–906. [PubMed]

4. Ortega F, Acerenza L, Westerhoff HV, Mas F, Cascante M. Product dependence and bifunctionality compromise the ultrasensitivity of signal transduction cascades. Proc Natl Acad Sci U S A. 2002;99:1170–1175. [PubMed]

5. Thattai M, van Oudenaarden A. Attenuation of noise in ultrasensitive signaling cascades. Biophys J. 2002;82:2943–2950. [PubMed]

6. Sauro HM, Kholodenko BN. Quantitative analysis of signaling networks. Prog Biophys Mol Biol. 2004;86:5–43. [PubMed]

7. Hornberg JJ, Bruggeman FJ, Westerhoff HV, Lankelma J. Cancer: A systems biology disease. Biosystems. 2006;83:81–90. [PubMed]

8. Huang CY, Ferrell JE., Jr Ultrasensitivity in the mitogen-activated protein kinase cascade. Proc Natl Acad Sci U S A. 1996;93:10078–10083. [PubMed]

9. Kholodenko BN. Negative feedback and ultra-sensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur J Biochem. 2000;267:1583–1588. [PubMed]

10. Blüthgen N, Herzel H. How robust are switches in intracellular signaling cascades? J Theor Biol. 2003;225:293–300. [PubMed]

11. Goldbeter A, Koshland DE., Jr An amplified sensitivity arising from covalent modification in biological systems. Proc Natl Acad Sci U S A. 1981;78:6840–6844. [PubMed]

12. Goldbeter A. A minimal cascade model for the mitotic oscillator involving cyclin and cdc2 kinase. Proc Natl Acad Sci U S A. 1991;88:9107–9111. [PubMed]

13. Igoshin OA, Goldbeter A, Kaiser D, Oster G. A biochemical oscillator explains several aspects of *Myxococcu*s *xanthu*s behavior during development. Proc Natl Acad Sci U S A. 2004;101:15760–15765. [PubMed]

14. Angeli D, Ferrell JE, Jr, Sontag ED. Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems. Proc Natl Acad Sci U S A. 2004;101:1822–1827. [PubMed]

15. Wolkenhauer O, Ullah M, Wellstead P, Cho K-H. A systems-and signal-oriented approach to intracellular dynamics. Biochem Soc Trans. 2005;33:507–515. [PubMed]

16. Csiksz-Nagy A, Kapuy O, Gyorffy B, Tyson JJ, Novk B. Modeling the septation initiation network (SIN) in fission yeast cells. Curr Genet. 2007;51:245–255. [PubMed]

17. Heinrich R, Neel BG, Rapoport TA. Mathematical models of protein kinase signal transduction. Mol Cell. 2002;9:957–970. [PubMed]

18. Nakabayashi J, Sasaki A. Optimal phosphorylation step number of intracellular signal-transduction pathway. Theor Biol. 2005;233:413–421. [PubMed]

19. Marhl M, Grubelnik V. Role of cascades in converting oscillatory signals into stationary step-like responses. Biosystems. 2007;87:58–67. [PubMed]

20. Millat T, Bullinger E, Rohwer J, Wolkenhauer O. Approximations and their consequences for dynamic modelling of signal transduction pathways. Math Biosci. 2006;207:40–57. [PubMed]

21. Legewie S, Schoeberl B, Blüthgen N, Herzel H. Competing docking interactions can bring about bistability in the MAPK cascade. Biophys J. 2007;93:2279–2288. [PubMed]

22. Cleland W W. The kinetics of enzyme-catalyzed reactions with two or more substrates or products II. Inhibition: Nomenclature and theory. Biochim. Biophys. Acta. 1963;67:173–187. [PubMed]

23. Salazar C, Höfer T. Kinetic models of phosphorylation cycles: a systematic approach using the rapid-equilibrium approximation for protein-protein interactions. Biosystems. 2006;83:195–206. [PubMed]

24. Gouzé JL. Positive and Negative circuits in dynamical systems. J Biol Systems. 1998;6:11–15.

25. Bhalla US, Iyengar R. Emergent properties of networks of biological signaling pathways. Science. 1999;283:381–387. [PubMed]

26. Levchenko A, Bruck J, Sternberg PW. Scaffold proteins may biphasically affect the levels of mitogen-activated protein kinase signaling and reduce its threshold properties. Proc Natl Acad Sci U S A. 2000;97:5818–5823. [PubMed]

27. Ferrell JE., Jr Tripping the switch fantastic: how a protein kinase cascade can convert graded inputs into switch-like outputs. Trends Biochem Sci. 1996;21:460–466. [PubMed]

28. Kholodenko BN, Kiyatkin A, Bruggeman FJ, Sontag E, Westerhoff HV, et al. Untangling the wires: a strategy to trace functional interactions in signaling and gene networks. Proc Natl Acad Sci U S A. 2002;99:12841–12846. [PubMed]

29. Kholodenko BN, Hoek JB, Westerhoff HV, Brown GC. Quantification of information transfer via cellular signal transduction pathways. FEBS Lett. 1997;414:430–434. [PubMed]

30. Santos SD, Verveer PJ, Bastiaens PI. Growth factor-induced MAPK network topology shapes Erk response determining PC-12 cell fate. Nat Cell Biol. 2007;9:324–330. [PubMed]

31. Kholodenko BN. Untangling the signalling wires. Nat Cell Biol. 2007;9:247–249. [PMC free article] [PubMed]

32. Gonze D, Goldbeter A. A model for a network of phosphorylation+dephosphorylation cycles displaying the dynamics of dominoes and clocks. J Theor Biol. 2000;210:167–186. [PubMed]

33. Murray AW, Kirshner MW. Dominoes and clocks: the union of two views of the cell cycle. Science. 1989;246:614–621. [PubMed]

34. Borghans JA, De Boer RJ, Segel LA. Extending the quasi-steady state approximation by changing variables. Bull Math Biol. 1996;58:43–63. [PubMed]

35. Ciliberto A, Capuani F, Tyson JJ. Modeling networks of coupled enzymatic reactions using the total quasi-steatdy state approximation. PloS Comput Biol. 2007;3:463–472.

36. Qiao L, Nachbar RB, Kevrekidis IG, Shvartsman SY. Bistability and oscillations in the Huang-Ferrell model of MAPK signaling. PLoS Comput Biol. 2007;3:1819–1826. [PMC free article] [PubMed]

37. Markevich NI, Hoek JB, Kholodenko BN. Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J Cell Biol. 2004;164:353–359. [PMC free article] [PubMed]

Articles from PLoS Computational Biology are provided here courtesy of **Public Library of Science**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |