Ultrasensitivity in signaling cascades revisited: Linking local and global ultrasensitivity estimations

Ultrasensitive response motifs, capable of converting graded stimuli into binary responses, are well-conserved in signal transduction networks. Although it has been shown that a cascade arrangement of multiple ultrasensitive modules can enhance the system’s ultrasensitivity, how a given combination of layers affects a cascade’s ultrasensitivity remains an open question for the general case. Here, we introduce a methodology that allows us to determine the presence of sequestration effects and to quantify the relative contribution of each module to the overall cascade’s ultrasensitivity. The proposed analysis framework provides a natural link between global and local ultrasensitivity descriptors and it is particularly well-suited to characterize and understand mathematical models used to study real biological systems. As a case study, we have considered three mathematical models introduced by O’Shaughnessy et al. to study a tunable synthetic MAPK cascade, and we show how our methodology can help modelers better understand alternative models.


Introduction
Sigmoidal input-output response modules are very well-conserved in cell signaling networks that might be used to implement binary responses, a key element in cellular decision processes.Additionally, sigmoidal modules might be part of more complex structures, where they can provide the nonlinearities which are needed in a broad spectrum of biological processes [1,2], such as multistability [3,4], adaptation [5], and oscillations [6].There are several molecular mechanisms that are able to produce sigmoidal responses such as inhibition by a titration process [7,8], zero-order ultrasensitivity in covalent cycles [9,10], and multistep activation processes -like multisite phosphorylation [11][12][13] or ligand binding to multimeric receptors [14].
Sigmoidal curves are characterized by a sharp transition from low to high output following a slight change of the input.The steepness of this transition is called ultrasensitivity [15].In general, the following operational definition of the Hill coefficient may be used to calculate the overall ultrasensitivity of sigmoidal modules: where EC10 and EC90 are the signal values needed to produce an output of 10% and 90% of the maximal response.The Hill coefficient n H quantifies the steepness of a function relative to the hyperbolic response function which is defined as not ultrasensitive and has n H = 1 (i.e. an 81-fold increase in the input signal is required to change the output level from 10% to 90% of its maximal value).Functions with n H > 1 need a smaller input fold increase to produce such output change, and are thus called ultrasensitive functions.
Global sensitivity measures such the one described by eq. 1 do not fully characterize s-shaped curves, y(x), because they average out local characteristics of the analyzed response functions.Instead, these local features are well captured by the logarithmic gain or response coefficient measure [16] defined as: Equation 2 provides local ultrasensitivity estimates given by the local polynomial order of the response function.

R(x) =
x y Equation 2 provides local ultrasensitivity estimates given by the local polynomial order of the response function.

MAP kinase cascades
Mitogen activated protein (MAP) kinase cascades are a well-conserved motif.They can be found in a broad variety of cell fate decision systems involving processes such as proliferation, differentiation, survival, development, stress response and apoptosis [17].They are composed of a chain of three kinases which sequentially activate one another, through single or multiple phosphorylation events.A thoughtful experimental and mathematical study of this kind of systems was performed by Ferrell and collaborators, who analyzed the steady-state response of a MAPK cascade that operates during the maturation process in Xenopus oocytes [18].They developed a biochemical model to study the ultrasensitivity displayed along the cascade levels and reported that the combination of the different ultrasensitive layers in a multilayer structure produced an enhancement of the overall system's global ultrasensitivity [18].In the same line, Brown et al. [19] showed that if the dose-response curve, F(x), of a cascade could be described as the mathematical composition of functions, fisi, that described the behavior of each layer in isolation (i.e, )) then the local ultrasensitivity of the different layers combines multiplicatively: R( In connection with this result, Ferrell showed for Hill-type modules of the form where the parameter EC50 corresponds to the value of input that elicits half-maximal output, and nH is the Hill coefficient), that the overall cascade global ultrasensitivity had to be less than or equal to the product of the global ultrasensitivity estimations of each cascade's layer, i.e n H ≤ n H,1 .nH,2 [20].
Hill functions of the form given by eq. 3 are normally used as empirical approximations of sigmoidal dose-response curves, even in the absence of any mechanistic foundation [2].However, it is worth noting that for different and more specific sigmoidal transfer functions, qualitatively different results could have been obtained.In particular, a supramultiplicative behavior (the ultrasensitivity of the combination of layers is higher than the product of individual ultrasensitivities) might be observed for left-ultrasensitive response functions, i.e. functions that are steeper to the left of the EC50 than to the right [21]).In this case, the boost in the ultrasensitivity emerges from a better exploitation of the ultrasensitivity "stored" in the asymmetries of the dose-response functional form (see [21] for details).
As modules are embedded in bigger networks, constraints in the range of inputs that each module would receive (as well as in the range of outputs that the network would be able to transmit) could arise.We formalized this idea in a recent publication introducing the notion of dynamic range constraint of a module's dose-response function.The later concept is a feature inherently linked to the coupling of modules in a multilayer structure, and resulted a significant element to explain the overall ultrasensitivity displayed by a cascade [21].Besides dynamic range constraint effects sequestration -i.e., the reduction in free active enzyme due to its accumulation in complex with its substrate-is another relevant process inherent to cascading that could reduce the cascade's ultrasensitivity [22][23][24].Moreover, sequestration may alter the qualitative features of any well-characterized module when integrated with upstream and downstream components, thereby limiting the validity of module-based descriptions [25][26][27].
All these considerations expose the relevance of studying the behavior of modular processing units embedded in their physiological context.Although there has been significant progress in the understanding of kinase cascades, how the combination of layers affects the cascade's ultrasensitivity remains an open-ended question for the general case.Sequestration and dynamic range constraints not only contribute with their individual complexity, but also usually occur together, thus making it more difficult to identify their individual effective contribution to the system's overall ultrasensitivity.
In the present work, we have developed a method to describe the overall ultrasensitivity of a molecular cascade in terms of the effective contribution of each module.In addition, said method allows us to disentangle the effects of sequestration and dynamic range constraints.We used our approach to analyze a recently presented synthetic MAPK cascade experimentally engineered by O'Shaughnessy et al. [28].
Using a synthetic biology approach O'Shaughnessy et al.
[28] constructed an isolated mammalian MAPK cascade (a Raf-MEK-ERK system) in yeast and analyzed its information processing capabilities under different rather well-controlled environmental conditions.They made use of a mechanistic mathematical description to account for their experimental observations.Their model was very similar in spirit to Huang-Ferrell's with two important differences: a) no phosphatases were included, and b) the creation and degradation of all species was explicitly taken into account.Interestingly, they reported that the multilayer structure of the analyzed cascades can accumulate ultrasensitivity supra-multiplicatively, and suggested that cascading itself and not any other process (such as multi-step phosphorylation, or zero-order ultrasensitivity) was at the origin of the observed ultrasensitivity.They called this mechanism, de-novo ultrasensitivity generation.As we found the proposed mechanism a rather appealing and unexpected way of ultrasensitivity generation, we wanted to further characterize it within our analysis framework.In particular, we reasoned that the methodology and concepts introduced in the present contribution were particularly well-suited to understand the mechanisms laying behind the ultrasensitivity behavior displayed by O'Shaughnessy cascade model.
The paper was organized as follows.First, we presented a formal connection between local and global descriptors of a module's ultrasensitivity for the case of a cascade system composed of N units.We then introduced the notion of Hill input's working range in order to analyze the contribution to the overall system's ultrasensitivity of a module embedded in a cascade.Next, we presented the O'Shaughnessy cascade analysis, in order to show the insights that might be gained using the introduced concepts and analysis methodologies.We concluded by presenting a summarizing discussion after which conclusions were drawn.

Linking local and global ultrasensitivity estimations
The concept of ultrasensitivity describes a module's ability to amplify small changes in input values into larger changes in output values.It is customary to quantify and characterize the extent of the amplification both globally, using the Hill coefficient n H defined in equation 1, and locally, using the response coefficient, R(I), as a function of the module's input signal I (equation 2), We found a simple relationship between both descriptions considering the logarithmic amplification coefficient A f a,b , defined as: A f a,b describes (in a logarithmic scale) the change produced in the output when the input varies from a to b values.For instance, A f a,b = 0.5 for an hyperbolic function evaluated between the inputs that resulted in 90% and 10% of the maximal output.In this case, the two considered input levels delimited the input range that should be considered for the estimation of the respective Hill coefficient n H ., We called this input interval: the Hill input's working range (see Fig 1A-B).-E).The X10 i and X90 i are the input values that take the "i" module when the last module (the second one) reach the %10 and %90 of it maximal Output (O 2,max ).when O 1,max >> EC50 2 (D) O 2,max equals the maximum output of module 2 in isolation, thus, X10 2 and X90 2 match the EC10 and EC90 of module 2 in isolation.Also the Hill working range of module 1 is located in the input region below EC501.On the other hand, when O 1,max EC50 2 (E) O 2,max is less than the maximum output of module 2 in isolation, thus, X10 2 and X90 2 differ from the EC10 and EC90 of module 2 in isolation.In this case the module's-1 Hill working range tends to be centered in values higher than its EC50, this will depends on module's-2 ultrasensitivity (see supplementary materials A) Taking into account eq. ( 4), the parameter n H could be rewritten as follows, Consequently, the Hill coefficient could be interpreted as the ratio of the logarithmic amplification coefficients of the function of interest and an hyperbolic function, evaluated in the corresponding Hill input's working range.
It is worth noting that the logarithmic amplification coefficient that appeared in equation 5 equaled the slope of the line that passed through the points (EC10, f (EC10)) and (EC90, f (EC90)) in a log-log scale.Thus, it was equal to the average response coefficient calculated over the interval [EC10, EC90] in logarithmic scale (see Fig 1c).Therefore, In particular, it can be appreciated that the module's Hill coefficient was related to the average response coefficient over the module's Hill input's working range, in units of a reference hyperbolic curve.

Ultrasensitivity in function composition.
We generalized the last result to cast the overall global ultrasensitivity level of a multitier cascade in terms of logarithmic amplification coefficients.We proceeded by first considering two coupled ultrasensitive modules, disregarding effects of sequestration of molecular components between layers.In this case, the expression for the system's dose-response curve, F , results from the mathematical composition of the functions, f i , which describe the input/output relationship of isolated modules i = 1, 2: Using eq. 5, the system's Hill coefficient n H resulted: where X10 i and X90 i delimited the Hill input's working range of the composite system, i.e. the input values for the i-layer so that the last layer (corresponding to i=2 in this case) reached the 10% and 90% of it maximal output level, respectively (see Fig 1D -1E).
It followed from equation 8 that the system's Hill coefficient nH could be written as the product of two factors, ν 1 and ν 2 , which characterized local average sensitivities over the relevant input region for each layer: The factor ν 2 in equation 8 was formally equivalent to the Hill coefficient of layer-2 but, importantly, now it was calculated using layer-1 Hill input's working range limits, X101 and X901, instead of the Hill working range limits of layer-1 in isolation, EC10 1 and EC90 1 .On the other hand, ν 1 was the amplification factor In this context, we dubbed the ν i coefficient: effective ultrasensitivity coefficient of layer-i, as it was associated to the effective contribution of layer-i to the system's overall ultrasensitivity.
For the more general case of a cascade of N modules we found that: This last equation, which hold exactly, showed a very general result.For the general case, the overall n H of a cascade could be understood as a multiplicative combination of the ν i of each module.The connection between the global and local ultrasensitivity descriptors, provided by equation 9, proved to be a useful tool to analyze ultrasensitivity in cascades, as it allowed assessing the effective contribution of each module to the system's overall ultrasensitivity.
According to this equation, Hill input's working range designated regions of inputs over which the mean local-ultrasensitivity value was calculated for each cascade level in order to set the system's n H .It was thus a significant parameter to characterize the overall ultrasensitivity of multilayer structures.In figures 1D-1E we illustrated, for the case of two composed Hill functions, to what extent the actual location of these relevant intervals depended on the way in which the cascade layers were coupled.The ratio between O 1,max and EC50 2 played a key role at this respect (see Fig 1).We showed how this parameter sets the Hill working ranges for the case of modules presenting different dose-response curves in Sup Mat. A. Importantly this analysis highlighted the impact of the detailed functional form of a module's response curve on the overall system's ultrasensitivity in cascade architectures.Local sensitivity features of the involved transfer functions were of the uttermost importance in this kind of setting and could be at the core of non-trivial phenomenology.For example a dose-response module with larger local ultrasensitivities than its overall global value might contribute with more ultrasensitivity to the system than the function's own n H (see Sup.Mat.A).

A 3-step analysis of cascade ultrasensitivity
In order to disentangle the different factors contributing to a cascade overall ultrasensitivity, we simultaneously considered three approximations of the system under study (see Fig 2).In a first step (Fig 2A ) we numerically computed the transfer function, f is , of each module in isolation and calculated the respective Hill coefficients n is .
We then studied the mathematical composition, F non−seq , of the isolated response functions (Fig 2B ): F non−seq represented the transfer function of the kinase cascade when sequestration effects were completely neglected.Following equation 9 effective ultrasensitivities, ν non−seq could be estimated and compared against the global ultrasensitivity that each module displayed in isolation, .Thus, this second step aimed to specifically analyze to what extent  In our analysis, we defined the output of a module and the input to the next one as the total active form of a species, including complexes with the next layer's substrates.However, we excluded complexes formed by same layer components (such as a complex between the phosphorylated kinase and its phosphatase), since these species are "internal" to each module.By doing this, we are able to consistently identify layers with modules (the same input/output definition was used by Ventura et al [25]).
In order to study the contribution of each layer to the system's ultrasensitivity, we proceeded to numerically compute the transfer function of each module in isolation and then calculate their respective n is H (column 2 in Table 1).The analysis of the mathematical composition of the three isolated transfer functions (step-2) allowed us to assess for the effects of module cascading.The explicit module coupling induced the existence of non-trivial Hill's input working ranges that were not captured in step-1 estimation.In this case, we found that the Fnon-seq Hill coefficient  resulted in n non−seq H = 3.91 (Table 1, column 3).This value was lower than the product of each module's Hill coefficient obtained in step 1 (n is 1 n is 2 n is 3 = 5.02).Thus, making the cascade to have a sub-multiplicative behavior.This result implied that inside the cascade arrangement at least one tier contributed less ultrasensitivity to the system than when considered acting in isolation.In Table 1 we displayed the values of the Hill coefficient for each layer when acting in isolation (n is ), and their effective ultrasensitivity coefficients (ν non−seq ) in columns 2 and 3 respectively.It can be seen that the resetting of the respective Hill working ranges produced a 13% and 11% decrease respectively in the effective contribution of the MAPKKK and MAPKK tiers to the overall system's ultrasensitivity.
Interestingly, we obtained that sequestration is not affecting the ultrasensitivity of this system, given that ν non−seq = ν seq for each cascade tier (columns 3 and 4 of Table 1).Moreover, we found that sequestration effects were actually negligible for the MAPKK and MAPK layers (see Supplementary materials B).
In order to probe the origin of the ultrasensitivity observed in the original cascade .The use of these models aimed to control for multi-step phosphorylation and complex formation effects on the resulting overall cascade ultrasensitivity, respectively.We considered rather illustrative to revisit these alternative models within our analysis framework.  .Because the cascade is not subject to multiple activation processes, competitive inhibition or zero-order ultrasensitivity (due to the absence of phosphatases), they claim that there is no other ultrasensitive source than the kinase-cascading architecture itself.The effective ultrasensitivity coefficients for this system are reported in Table 2.It can be seen that sequestration is not an important factor in this system.Moreover, the reported ν non−seq values suggest that the system's ultrasensitivity is mainly provided by ERK ultrasensitivity.Synthesis and degradation happened to be the key factors to understand the origin of the ERK layer's ultrasensitivity.This layer (scheme in Fig 4A ) was in fact mathematically analogous to a covalent cycle (scheme in Fig 4B ) because there was an implicit channel from the activated protein towards its inactive form via the degradation of the active protein and the production of the inactive form.Given that degradation was a linear reaction with respect to the amount of activated protein, its mathematical description was equivalent to a dephosphorylation reaction operating in a first order regime (Fig 4B ).Thus, the one-step system depicted in Fig 4 could in fact be described by a Goldbeter-Koshland (G-K) [9] function with We plotted in Fig 4C the steady state transfer function of the ERK module in isolation and the corresponding centered G-K function (see Supplementary Materials C).A clear agreement between both functions can be appreciated.Both curves present n H = 1.76, which is the ultrasensitivity contributed to the entire cascade.In the light of these results we conclude that the single-step cascade's ultrasensitivity did not come from a cascading effect but from a 'hidden' first-order ultrasensitivity process in the ERK layer.  .This system presented an overall Hill coefficient of 2.7, which is less than that value found for the original system (n H = 3.91).
In order to identify the origin of this decrease in ultrasensitivity, we calculated the corresponding effective ultrasensitive coefficients shown in Table 3.
We verified that the Hill function is not contributing as much ultrasensitivity as the original system's Raf-MEK module.The reason is that even the dose-response of active MEK and the Hill function appear to be similar, there are strong dissimilarities in their local ultrasensitivity behavior (see Fig 5B).This is particularly true for low input values, where the Hill's input working range is located.In this region, the active MEK curve presents local ultrasensitivity values larger than the Hill function counterpart, thus the replacement by a hill function produce a reduction in the Hill coefficient In this way, despite the high-quality of the fitting adjustment Residual Standard Error=2.6), the Hill function approximation introduced non-trivial alterations in the system's ultrasensitivity as a technical glitch.

Discussion
There have been early efforts to interpret cascade system-level ultrasensitivity out of the sigmoidal character of their constituent modular components.Usually they have focused either on a local or a global characterization of ultrasensitivity features.For instance, Brown et al. [19] had shown that the local system's ultrasensitivity in cascades equals the product of the local ultrasensitivity of each layer.In turn, from a global ultrasensitivity perspective, Ferrell [20] pointed out that in the composition of two Hill functions, the Hill coefficient results equal or less than the product of the Hill coefficient of both curves (n H ≤ n H,1 n H,2 ).
In this contribution we have found a mathematical expression (Equation 6) that linked both, the local and global ultrasensitivity descriptors in a fairly simple way.Moreover we could provide a generalized result to handle the case of a linear arrangement of an arbitrary number of such modules (Equation 9).Noticeably, within the proposed analysis framework, we could decompose the overall global ultrasensitivity in terms of a product of single layer effective ultrasensitivities.These new parameters were calculated as local-ultrasensitivity values averaged over meaningful working ranges (dubbed Hill's input working ranges), and permitted to assess the effective contribution of each module to the system's overall ultrasensitivity.Of course, the reason why we could state an exact general equation for a system-level feature in terms of individual modular information was that in fact system-level information was used in the definition of the Hill's input working ranges that entered Equation 9.The specific coupling between ultrasensitive curves, set the corresponding Hill's input working ranges, thus determining the effective contribution of each module to the cascade's ultrasensitivity.This process, which we called Hill's input working range setting, has already been noticed by several authors [20,29,30,23], but as far as we know this was the first time that a mathematical framework, like the one we present here, has been proposed for it.
The value of the obtained expression (Equation 9) resides in the fact that not only it captured previous results, like Ferrell's inequality, but also that it threw light about the mechanisms involved in the ultrasensitivity generation.For instance, the existence of supramultiplicative behavior in signaling cascades have been reported by several authors [23, 28] but in many cases the ultimate origin of supramultiplicativity remained elusive.Our framework naturally suggested a general scenario where supramultiplicative behavior could take place.This could occur when, for a given module, the corresponding Hill's input working range was located in an input region with local ultrasensitivities higher than the global ultrasensitivity of the respective dose-response curve.
In order to study how multiple ultrasensitive modules combined to produce an enhancement of the system's ultrasensitivity, we have developed an analysis methodology that allowed us to quantify the effective contribution of each module to the cascade's ultrasensitivity and to determine the impact of sequestration effects in the system's ultrasensitivity.This method was particularly well suited to study the ultrasensitivity in MAP kinase cascades.
We used our methodology to revisit O'Shaughnessy et al. tunable synthetic MAPK system [28] in which they claim to have found a new source of ultrasensitivity called: de novo ultrasensitivity generation.They explained this new effect in terms of the presence of intermediate elements in the kinase-cascade architecture.
We started analyzing the MAPK cascade.We found that sequestration was not affecting the system's ultrasensitivity and that the overall sub-multiplicative behavior was only due to a re-setting of the Hill input's working range for the first and second levels of the cascade.
Then, to investigate the origin of the claimed de novo ultrasensitivity generation mechanism, we applied our framework on the single-step phosphorylation cascade.We found that the system's ultrasensitivity in the single-step cascade came only from the contribution of the last module, which behaved as a Goldbeter-Koshland unit with kinases working in saturation and phosphatases in a linear regime.Therefore the ultrasensitivity in its single-step cascade was not generated by the cascading itself, but by the third layer, which itself was actually ultrasensitive.
Finally we analyzed the auxiliary model considered by O'Shaughnessy et al. in which the Raf and MEK layers were replaced by a Hill function that is coupled to the ERK layer.In this case, ,even the original Estradiol-MEK input-output response curve could be fairly well fitted and global ultrasensitivity features were rather well captured, the replacement by a Hill function produce a strong decrease in the systems ultrasensitivity.We found that the functional form of the Hill function failed to reproduce original local ultrasensitivity features that were in fact the ones that, due to the particular Hill working range setting acting in this case, were responsible for the overall systems ultrasensitivity behavior.The analyzed case was particularly relevant, as provided an illustrative example that warned against possible technical glitches that could arise as a consequence of the inclusion of approximating functions in MAPK models.

Conclusion
The study of signal transmission and information processing inside the cell has been, and still is, an active field of research.In particular, the analysis of cascades of sigmoidal modules has received a lot of attention as they are well-conserved motifs that can be found in many cell fate decision systems.In the present contribution we focused on the analysis of the ultrasensitive character of this kind of molecular systems.
We presented a mathematical link between global and local characterizations of the ultrasensitive nature of a sigmoidal unit and generalized this result to handle the case of a linear arrangement of such modules.In this way, the overall system ultrasensitivity could define in terms of the effective contribution of each cascade tier.Based on our finding, we proposed a methodological procedure to analyzed cascade modular systems, in particular MAPK cascades.
We used our methodology to revisit O'Shaughnessy et al. tunable synthetic MAPK system [28].In which they claim to find a new source of ultrasensitivity called: de novo ultrasensitivity generation, which they explained in terms of the presence of intermediate elements in the kinase-cascade architecture.With our framework we found that the ultrasensitivity did not come from a cascading effect but from a 'hidden' first-order ultrasensitivity process in the one of the cascade's layer.
From a general perspective, our framework serves to understand the origin of ultrasensitivity in multilayer structures, which could be a powerful tool in the designing of synthetic systems.In particular, in ultrasensitive module designing, our method can be used to guide the tuning of both the module itself and the coupling with the system, in order to set the working range in the region of maximal local ultrasensitivity.A The effect of the Hill's input working range in multitiered systems.
The Hill's input working range delimits the region of inputs over which the mean localultrasensitivity value is calculated (equation 9).It is thus a significant parameter to get insights about the overall ultrasensitivity of multilayer structures.In what follows, we show that the actual location of this relevant interval depends on the way in which cascade layers are coupled.Let's start by considering two coupled ultrasensitive modules.Two different regimes could be identified depending whether the upstream module's maximum output was or wasn't large enough to fully activate the downstream unit: a) In the first case i.e. when O 1,max EC50 2 (see Fig 1D), X10 2 and X90 2 are equal to the EC10 2 and EC90 2 levels respectively.Therefore, when coupled to module-1, the Hill input's working range of module-2 would not differ from the isolated case, and would equal the Hill coefficient of this module acting in isolation: i.e ν 2 = n is 2 .In addition, it can be seen that the Hill input's working range of module-1 tends to be located at the low input-values region for increasing levels of the ratio O 1,max /EC50 2 .In this region the response coefficient of the Hill functions achieve the highest values, with R 1 ≈ n is 1 (see Fig 1C), thus, when calculating the average logarithmic gain, module's-1 working range in a region with low local ultrasensitivity.Although we show that the submultiplicativity occurs in the limit of O 1,max EC50 2 , the same argument is still valid for O 1,max EC50 2 .
Of course, the ultimate consequences in the coupling of two ultrasensitive modules will depend on the particular mathematical details of the transfer functions under consideration.In this way, a completely qualitatively different behavior could be found for a system composed of two modules characterized by Golbeter-Koshland, GK, response functions [2].GK functions appear in the mathematical characterization of covalent modification cycles (such as phosphorylation-dephosphorylation), ubiquitous in cell signaling, operating in saturation.For cases where the phosphatases, but not the kinases, work in saturation, GK functions present input regions with response coefficients higher than their overall n H (see Fig 2A-C).Their Hill input's working range are thus located in the region of greatest local ultrasensitivity, these functions are able to contribute with more effective ultrasensitivity than their global ultrasensitivity.Therefore, cascades involving GK functions may exhibit supra-multiplicative behavior.For this kind of systems, Fig 2D shows that, under regime (a) (when the module's-2 EC50 is much lower than the GK maximal output level, O1 max) the module's-1 Hill input's working range is set in its linear regime (R 1 = 1), and the GK function does not contribute to the overall system's ultrasensitivity.On the other hand, Fig 2E shows that the O 1,max to EC50 2 relation can be tuned in order to set module's-1 Hill input's working range in its most ultrasensitive region, producing an effective ultrasensitivity contribution, ν 2 , even larger than the ultrasensitivity of the GK curve in isolation (i.e.ν 2 ≥ n 2 ), resulting in supra-multiplicative behavior n H = ν 1 .ν 2 > n is 1 .n is 2 .Our analysis highlights the impact of the detailed functional form of a module's response curve on the overall system's ultrasensitivity in cascade architectures.Local sensitivity features of the involved transfer functions are of the uttermost importance in this kind of setting and could be at the core of non-trivial phenomenology.
havior was only due to a resetting of the Hill input's working range for the first and second levels of the cascade.Transfer function of each isolated layers (Is), compared with the transfer function that this module is actually sustaining in the cascade (Seq) and in the composition of isolated functions (Non-Seq), (A-C).Respective response coefficient curves (D-F).The blue dashed vertical lines show the X10 i and X90 i of each layer of the original cascade, while the red solid vertical lines show the X10 i and X90 i of each layers in the system given by the composition of the modules in isolation (F non−seq , see eq [6]).It worth noting that the response curves that each module sustain in the non-sequestration scenario (panels B-C) will coincide with the isolated curves, with the exception that are limited in the spanned input region.

C Goldbeter-Koshland function
The Goldbeter-Koshland function [1] is defined as In order to center the G-K function, we multiply the independent variable for a scale factor α, G(α, x, K 1 , K 2 ), where α was set in order to make the EC50 of G-K function coincides with the EC50 of ERKpp curve

Fig. 1 :
Fig. 1: Hill function dose-response, and its composition.Schematic representation of Hill type dose-response curve in (A) log-linear scale and (B) loglog scale, and (C) respective local ultrasensitivity.The EC10 and EC90 are the signal values needed to produce an output level of 10% and 90% of the maximal response (O max ), and the "Hill working range" is the input range relevant for the calculation of the system's n H ill. In Hill functions, inputs values much smaller than its EC50 produce local sensitivities around its Hill coefficient.Schematic response function diagrams for the composition of two Hill type ultrasensitive modules (D-E).The X10 i and X90 i are the input values that take the "i" module when the last module (the second one) reach the %10 and %90 of it maximal Output (O 2,max ).when O 1,max >> EC50 2 (D) O 2,max equals the maximum output of module 2 in isolation, thus, X10 2 and X90 2 match the EC10 and EC90 of module 2 in isolation.Also the Hill working range of module 1 is located in the input region below EC501.On the other hand, when O 1,max EC50 2 (E) O 2,max is less than the maximum output of module 2 in isolation, thus, X10 2 and X90 2 differ from the EC10 and EC90 of module 2 in isolation.In this case the module's-1 Hill working range tends to be centered in values higher than its EC50, this will depends on module's-2 ultrasensitivity (see supplementary materials A) where X a,b denoted the mean value of the variable x over the range [a,b].This last equation explicitly linked the local and global ultrasensitivity descriptions.

Fig. 2 :
Fig. 2: Modular and system representation of a MAP kinase cascade.Each layer in isolation is composed by single or multiple covalent cycles, which dose-response curves can be ultrasensitive by zero-order mechanisms and/or multi-activation processes (A).The cascade transfer function, in a scenario in which sequestration is not taken into account (F non−seq ) can be obtain by the mathematical composition of each module's transfer functions f is i acting in isolation (B).When the sequestration effect is taken into account, the layers embedded in the MAP kinase cascade may have a different dose-response curve from the isolated case (C).

2. 4
Ultrasensitivity in O'Shaughnessy et al. model A sketch of the O'Shaughnessy et al.MAPK cascade is shown in Fig 3A.

Fig. 3 :
Fig. 3: O'Shaughnessy et al. cascades scheme.Dual-step cascade scheme (A), Single-step cascade scheme (B), and Dual-step cascade scheme with Raf-MEK system replaced by a Hill Function (C).In each case, Estradiol is the input and the most phosphorylated state of ERK is the output.

(
sketched in Fig 3A), O'Shaughnessy et al. analyzed a couple of alternative computational models.In an auxiliary model (a) they replaced dual-step by single-step phosphorylation layers (see Fig 4B).In a second auxiliary model (b) they aggregated the Raf and MEK levels of the cascade into a single Hill equation that mimicked the dose-response of the total active MEK in the original model (see Fig 4C)

Fig. 4 :
Fig. 4: Equivalence between a single-step layer in O'Shaughnessy model and a covalent cycle.O'Shaughnessy et al. single-Step layer scheme (A) and equivalent covalent cycle scheme (B).Steady state transfer function of ERK layer in isolation of the O'Shaughnessy single-step cascade in blue dashed line (C), compared to a centered Goldbeter-Function with equivalent parameters in red solid line (K 1 = 0.04 and K 2 = 1000, see Sup.Mat.C )

2. 4 . 2
Analysis of auxiliary model (b): the devil is in the details In Fig 3C, we sketched the system resulting from the composition of a Hill function with the ERK layer analyzed by O'Shaughnessy et al.The Hill parameters were worked out by fitting the active MEK dose-response by a Hill function (Fig 5A)

Fig. 5 :
Fig. 5: Fitting by a Hill function may neglect relevant behaviors.Dose-response curve of active MEK in O'Shaughnessy model compared with its fit by a Hill function (A).Respective response coefficient (B).It can be seen that as the dose-response of active MEK and the Hill function appear to be similar, there are strong dissimilarities in their local ultrasensitivity downstream constraints on a signaling module's ultrasensitivity.Physical biology.2014 Oct 14;11(6):066003.22. Bluthgen N, Bruggeman F J, Legewie S, Herzel H, Westerhoff H V and Kholodenko B N Effects of sequestration on signal transduction cascades.FEBS Journal 2006 273 895-906 23.Racz E and Slepchenko B M On sensitivity amplification in intracellular signaling cascades Phys.Biol.2008 5 036004-12 24.Wang G, Zhang M. Tunable ultrasensitivity: functional decoupling and biological insights.Scientific Reports.2016 6 20345 25.Ventura A C, Sepulchre J A and Merajver S D A hidden feedback in signaling cascades is revealed, PLoS Comput Biol.2008 4(3):e1000041 26.Del Vecchio D, Ninfa A and Sontag E Modular cell biology: retroactivity and insulation.Mol Sys Biol 2008 4 161 p161.27.Ventura A C, Jiang P, Van Wassenhove L, Del Vecchio D, Merajver S D, and Ninfa A J Signaling properties of a covalent modification cycle are altered by a downstream target PNAS 2010 107 (22) 10032-10037 28.O'Shaughnessy E C, Palani S, Collins J J, Sarkar C A Tunable signal processing in synthetic MAP kinase cascades Cell 2011 7 144(1):119-31 29.Bluthgen N and Herzel H MAP-Kinase-Cascade: Switch, Amplifier or Feedback Controller? 2nd Workshop on Computation of Biochemical Pathways and Genetic Networks -Berlin: Logos-Verlag 2001 55-62 30.Bluthgen N and Herzel H How robust are switches in intracellular signaling cascades Journal of theoretical Biology 2003 225 293-300

Fig. 7 :
Fig. 7: O'Shaughnessy model dose-response analysis.Transfer function of each isolated layers (Is), compared with the transfer function that this module is actually sustaining in the cascade (Seq) and in the composition of isolated functions (Non-Seq), (A-C).Respective response coefficient curves (D-F).The blue dashed vertical lines show the X10 i and X90 i of each layer of the original cascade, while the red solid vertical lines show the X10 i and X90 i of each layers in the system given by the composition of the modules in isolation (F non−seq , see eq[6]).It worth noting that the response curves that each module sustain in the non-sequestration scenario (panels B-C) will coincide with the isolated curves, with the exception that are limited in the spanned input region.

Table 1 :
3-step analysis of O'Shaughnessy et al. dual-step cascade.Column 1: Hill coefficients of active species dose-response curves (respect to estradiol).Column 2: Hill coefficients of the modules in.Column 3: Effective ultrasensitivities of the system given by the composition of the modules in isolation, which is represented by F non−seq (equation [6]) where no sequestration effects are present.Column 4: Effective ultrasensitivities in the original cascade.

Table 3 :
Ultrasensitivity in O'Shaughnessy et al. dual-step cascade when the Raf-ERK module is replaced by a Hill function.Column 1: Hill coefficients of active species dose-response curves (respect to estradiol).Column 2: Hill coefficients of each modules in isolation, which coincide with the dose-response curves in the original cascade as there is no sequestration effects.Column 3: Effective ultrasensitivity of each layers when they are coupled as a cascade.