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.
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.
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 ηi are ratios of total amounts of proteins. εi is the ratio of total phosphatase over total targeted protein. ηi is defined as the total targeted protein in one cycle over the corresponding amount in the next cycle in the cascade, or, equivalently, the ratio of total kinase over total targeted protein. The parameter µi is the ratio of the kinetic rates of product formation in both the activation and the inactivation reactions (see reactions in Equation 1). 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 xi = yi*+ci+1, which is the natural slow variable describing the total kinase i available at a given time for the phosphorylation in cycle i+1. This reduction is only valid if the total phosphatase in the cycle is much lower than the total targeted protein, i. e. , in the limit εi«1. The other parameters must satisfy µi ηi∼εi. The dynamics of xi is described by the differential equation: (6) with the following conservation equation from which yi has to be extracted: (7) x0 = S is the normalized input signal and yn+1 = 0. In Equation 6, Vi = (k′i µi ηi) / (ε k′) and V′i = (εi k′i) / (ε k′), where ε k′ is a typical number representing the set of εi k′i (i = 1, …, n), e. g. the arithmetical or the geometrical average over this set. In the conservation equation (Equation 7), the notation O (εi) is just a reminder that this equation is written in the lowest order in εi, as is also the case for the differential equation for xi. In Text S1, we discuss an improvement of this conservation equation which takes into account the first correction in εi. Although this extension does not alter the new properties discussed below, its numerical integration is easy and it increases the accuracy of the approximation. 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 yi+1, in contrast to the GK-like model. This function has the appearance of an effective Michaelis-Menten coefficient K′eff, i = K′i (1+yi+1/Ki+1), which is a typical way to indicate competitive inhibition in enzyme kinetics [22]. In the context of activation-inactivation cycles, a similar type of equation was obtained by Salazar and Höfer in the systematic study of a single cycle taking into account the competition between kinase and phosphatase to bind the same target protein [23]. In that case, an effective Michaelis-Menten coefficient appears also in the negative term of Equation 4, but with the form K′eff, i = (1+yi/Ki). In our study, however, the competition is induced by the next substrate yi+1, and this precisely describes a negative feedback from cycle i+1 on cycle i: the higher the level of xi+1, the smaller yi+1 and, therefore, the larger the value of the negative term in Equation 6. This modified denominator reflects the influence of the downstream step on the state variables of one given cycle. It is not a detail of the formalism. It has consequences upon the dynamics and on the properties of the signaling pathway, as it will be demonstrated in the following sections. Moreover, we will see that, since our system arises from a controlled approximation of the mechanistic model, the dynamics of both models can be made comparable. In the limit ηi∼εi«1, one retrieves the simple conservation law xi+yi≈1. However, we note that even in that limit and due to K′eff, i, our resulting system is not equivalent to the GK-like model. Notice that ηi∼εi«1 is the closest we can be to the hypothesis behind the GK-like model, where it is considered that the concentration of the targeted protein is in large excess over those of the converting enzymes. In our description, the converting enzymes for unit i are E′iT (phosphatase) and Yi−1, T (kinase). Taking the limit ηi«1 together with the fact that the targeted protein of each cycle is the activating protein of the next one, results in increasing protein concentrations as the cascade proceeds. Even though this is not the usual condition in signaling cascades, examples could arise where this limit is suitable. As a possible relevant example, the concentrations reported for the MAPK cascade go from nM in the first unit to µM in the second and third ones [8]. 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 ηi∼1, for all or for some index i, as long as µi ηi∼εi, which results in the limit µi∼εi«1. Since µi = ki/k′i, this limit requires that the phosphatase of that cycle be much more active that the corresponding kinase. In this limit however, the conservation law remains as expressed in Equation 7 and no further simplification can be made. As a result, in this limiting case, the first term in Equation 6 depends on the variable describing the previous step in a different (and more complicated) way compared to Equation 4. 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 µi∼εi½ and ηi∼εi½ for all or for some index i. 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 xi 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 yi*, the level of active protein, obtained from xi in Equations 6–7 (see Text S4 for a comparison between variables xi and yi*). In this section, each parameter in the reduced mechanistic model is considered to be homogeneous throughout the chain, i. e. , the parameters do not depend on the index i characterizing the position of a particular unit in the chain. The homogeneity assumption implies that Vi≡V = µη/ε and V′i≡V′ = 1. Parameter S indicates the level of input stimulation the chain receives. The parameter K′ is chosen by considering the relation K′ = K/µ. We have performed numerical simulations with other parameter relationships and the properties reported below are not critically dependent upon that choice. The control parameters are, then, V, K, µ, and η. Since V′ = 1, the range of V values of interest lies around 1. The initial condition for all the numerical simulations considered is, at t = 0, xi = 0 (and yi = 1) for every i. 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 K″ and take the values of 1,0. 25, and 0. 25, respectively. 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 E1T, the total amount of kinase for the first unit in the cascade (corresponding to MAPK kinase kinase in B). E1T, 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 y1*, y2*, and y3*, respectively). For completeness, we are also including the corresponding outcomes obtained by the GK-like model (dotted lines). In order to compare the steepness in the responses, we have computed the apparent Hill coefficient nH ([8]) for each curve, as indicated in the legend. As expected, nH increases through the chain. Moreover, nH is also considerably reduced when comparing GK-like model' s predictions with the predictions of both mechanistic and reduced mechanistic models (which are, as already mentioned, undistinguishable). As explained in the section dealing with stimulus-response curves, these differences could be due to the fact that both the mechanistic and the reduced mechanistic descriptions take into account “sequestration” in the enzymatic reactions [3]. 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 nH are not modified by removing the εi correction in the conservation law or, even by, removing the ηi correction as well (i. e. , using a conservation law of the form xi+yi = 1). These results strongly indicate the robustness of the new equations regarding the “ultrasensitivity” characteristics of the cascade. 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 y3* (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 y2* is actually a decrease due to a higher sequestration of free y2* by the next step in the chain caused, in turn, by the increased demand of y3. This result is exhibited by both cascades in Figure 7 (the one involving only single-phosphorylation and the one with double-phosphorylation in units 2 and 3) and we call it “reverse” stimulus-response curves. As stated before, this result is obtained with both the mechanistic and the reduced mechanistic descriptions, with realistic parameters associated with a well studied signaling pathway, such as MAPK. 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 y2* 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 y1* (black filled line), but this fact does not have to be generalized (see Text S6). The dotted horizontal lines in Figures 7A and 7B are the GK-like prediction for the response curve y2*: within that phenomenological description, a particular level in the cascade is not at all influenced by what happens in a downstream unit. However, this well known property of unidirectional influence in a signaling chain, which is embodied by the appellation of “cascade”, is shown here not to be guaranteed in general signaling cascades. 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/∂xi+1 indicates that the (i+1) th level of the cascade exerts a negative effect in what concerns variable xi. This effect (what we have called “negative feedback”) is intrinsic, as opposed to “explicit” negative feedback which is sometimes considered in models of signaling pathways [9], [12], [13]. MRA is, then, an appropriate approach to test this bidirectional structure and to estimate the relative strength of the backward interaction, as compared with the forward coupling in a signaling cascade. As a result of applying MRA, a matrix of local response coefficients r is obtained. An element rij in this matrix describes how the state of the variable associated with module j directly influences the state of the variable associated with module i. More precisely, a response coefficient rij lower/greater than 1 means that a relative change in module j is attenuated/amplified in module i by a factor rij (i. e. , Δxi/xi = rij Δxj/xj). A zero response coefficient indicates no direct effect between the involved modules, whereas a negative response coefficient means inhibition. In this way, the matrix r provides an interaction map to characterize the type and strength of the interactions between the modules in a cellular regulatory network. Indeed, if the rate of change of variable xi is denoted by the function fi, it easily can be shown that: (8) meaning that rij corresponds to a scaled version of the Jacobian matrix ∂fi/∂xj (evaluated in the steady state). Moreover, it was proven that the local response matrix r can be obtained from another matrix named global response matrix, Rp, that has the advantage of being accessible experimentally [28]. For example, the element (i, j) of this matrix can be obtained by perturbing a parameter pj affecting only module j and computing the relative changes induced on the steady state of xi, namely (Δxi/xi) /Δpj. For more details about the broad scope of the method, we refer the reader to the cited reference and references therein. 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 xi. Figure 8A contains the matrix of local response coefficients r. This matrix was obtained both by direct computation of the scaled Jacobian matrix (Equation 8) and by simulating experimental perturbations to the cascade, then computing the global response matrix, Rp, and finally obtaining r, as described previously [28] (details of second calculation not shown). Using MRA, the “theoretical” and “experimental” outputs were in perfect agreement and the results are displayed in Figure 8A. The structure of matrix r is tridiagonal, meaning that the first level in the cascade does not directly influence the third one (r31 = 0), and viceversa (r13 = 0). Coefficients r21 and r32 are positive, representing the positive effect of each level in the cascade to the subsequent one. Interestingly, r12 and r23 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 r21, r32, r12, and r23 versus the parameter E1T. r31 and r13 are zero (data not shown), r21 and r32 are positive, and r12 and r23 are negative, throughout the range where E1T was varied. Depending on the value of E1T, each of the nonzero rij 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. |r12/r21|, could exhibit large variations. Similar curves have been reported in the literature for signaling cascades [29], but lacking the information about r12 and r23, which have always been considered to be zero in previous papers. 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 xi (matrix r (xi) in Figure 8A) shows strong differences when compared to the one computed for the “free” enzyme variables yi* (matrix r (yi*) in Figure 8C. Although an explicit set of differential equations is not written for the variables yi*, the matrix r (yi*) can be calculated using the “experimental”' method described in the literature [28]. The result in Figure 8C is the average of four outputs and the corresponding error (standard error of the mean) is lower than 4%. As indicated in the reconstructed topology close to the matrix, r12 and r23 are now positive (as are r21 and r32), r31 is zero, and r13 is negative, indicating an inhibitory coupling from variable y3* to variable y1*. The matrix r (yi*) is consistent with the results in Figure 7A (and also those in Text S6): in other words, the response in y2* goes in the same direction as the one in y3* (whereas plotting variables xi indicates a decrease in x3 and an increase in x2, data not shown). Experimental data concerning the application of MRA to the MAPK cascade are now available in the literature [30], showing non zero r21 and r32 coefficients (and also non zero r31 and r13 coefficients). The interpretation of non zero r31 and r13 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 r12 and r23 coefficients was, at least regarding r23 and based on experimental evidence, that not only is y2* able to phosphorylate y3, but y3*can phosphorylate y2 as well [30], [31]. Our results however, suggest that the non vanishing backward coefficients (r12, r13, r23) can be accounted for, at least partly, by the natural “implicit” feedback which can exist in a signaling cascade. A quantitative correlation between these recent experimental results and our predictions is not possible at this time. In the published experiments, the MAPK cascade is not isolated but embedded in the complex cellular machinery; therefore, the measured connectivities could involve proteins external to the cascade itself and it would be premature to establish the connection with our simplified model. Nevertheless, the work in [30] suggests a direction for the type of experiments that could validate our results.
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 Yi* is not directly the kinase of the next reaction. Instead, we studied the case where Yi* activates that kinase. This scheme was suggested by the work of Goldbeter [12]. The resulting equation (see Text S7) is fundamentally different from the GK-like model. In reality, no set of assumptions can give rise to the GK-like model as a limiting case of our model. 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 ith cycle plus the amount of this protein which is captured by the next inter-converting cycle. The idea of working with a mixed variable xi can be further generalized by considering the “total” variable corresponding to the total amount of activated enzyme found not only as free molecules or bound to the next substrate, but also complexed with the reverse enzyme E′i. In fact, this choice is the key ingredient of the method called the “total” quasi-steady state approximation (tQSSA) which has been proved to be a simple but most efficient extension of the standard QSSA [34]. The application of this extended framework to the description of the signaling cascade of Figure 1 is concerned with our current research. In the same context, other authors have recently applied the tQSSA method to the study of small networks of GK cycles [35]. These systems do not form cascades, but involve a more complicated coupling between the units. Nevertheless, their results show that indeed the tQSSA method is successful in obtaining a reduced set of equations, with one variable per cycle, which faithfully reproduces the dynamics of the network for a large range of system parameters. 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 µi ηi∼εi. The study of the performance of the new approximation (Figure 2 and the corresponding computed errors) makes it clear that even when those conditions are satisfied only moderately, the new model is still robust in approximating the complete description. As an example, we have computed a 5% error for ε = 0. 1, η = 1, and µ = 0. 5 (meaning µη∼5ε). Moreover, we have observed that the steady state predictions of the reduced model are highly accurate. Therefore the properties of signaling cascades we are unveiling thanks to the new reduced model, are not restricted by a tight relationship between concentrations and reaction rates hard to achieve in in vivo or in vitro conditions. 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 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.
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.
computational biology/signaling networks
biochemistry/cell signaling and trafficking structures
biochemistry/theory and simulation
cell biology/cell signaling
A Hidden Feedback in Signaling Cascades Is Revealed