Retroactivity induced operating regime transition in an enzymatic futile cycle

Activated phosphorylation-dephosphorylation biochemical reaction cycles are a class of enzymatic futile cycles. A futile cycle such as a single MAPK cascade governed by two underlying enzymatic reactions permits Hyperbolic (H), Signal transducing (ST), Threshold-hyperbolic (TH) and Ultrasensitive (U) operating regimes that characterize input-output behaviour. Retroactive signalling caused by load due to sequestration of phosphorylated or unphosphorylated form of the substrate in a single enzymatic cascade without explicit feedback can introduce two-way communication, a feature not possible otherwise. We systematically characterize the operating regimes of a futile cycle subject to retroactivity in either of the substrate forms. We demonstrate that increasing retroactivity strength, which quantifies the downstream load, can trigger five possible regime transitions. Retroactivity strength is a reflection of the fraction of the substrate sequestered by its downstream target. Remarkably, the minimum required retroactivity strength to evidence any sequestration triggered regime transition demands 23% of the substrate bound to its downstream target. This minimum retroactivity strength corresponds to the transition of the dose-response curve from ST to H regime. We show that modulation of the saturation and unsaturation levels of the enzymatic reactions by retroactivity is the fundamental mechanism governing operating regime transition.


Introduction
Enzymatic cascades or futile cycles consisting of phosphorylation-dephosphorylation biochemical reaction cycles, are crucial, ubiquitously conserved, building-blocks of cellular signalling networks [1,2]. An enzymatic futile cycle employs phosphorylation and dephosphorylation reactions, respectively catalysed by kinase and phosphatase, to enable transition of a protein substrate between its two forms, namely inactive and active. Enzymatic cascades impart important properties like responsiveness, robustness, specificity onto a signalling response [3,4], weak signal amplification [5], signal speed acceleration [6], filter out noise in signal [7][8][9]. One such well-known enzymatic cascade is the Raf/MEK/ERK MAPK cascade, a key signal amplifier and a modulator of pro-survival and pro-apoptotic signalling pathways [10][11][12]. Aberrant functioning of this cascade has been implicated in many diseases such as cancer [13,14]. Detailed  MAPK cascade can therefore offer useful insights in designing therapeutic strategies for combating certain diseases. Activation behaviour of a futile cycle as a response to a stimulus of certain strength and their dynamic evolution have traditionally been characterized by systematically studying the dose-response curves permitted by the cascade [15][16][17]. Dose-response curve or the input-output characteristic of the cycle at steady-state is a map of the abundances of the input kinase and of the active protein (output) of the cascade [18]. Based on the qualitative nature of the dose-response curve, dictated by the saturated/unsaturated state of the two enzymatic reactions, the activation behaviour of futile cycle have been classified into four distinct operating regimes, viz., Hyperbolic (H), Signal transducing (ST), Threshold-hyperbolic (TH), Ultrasensitive (U), each of which display different signal processing capabilities [18,19]. Operating regimes of the MAPK cascades juxtaposed with patient-stratification data have recently been considered in disease prognostics [20]. Recently, a hybrid deterministic-stochastic approach constrained by experimental ensemble data was used for predicting and characterising the input-output behaviour of a single MAPK cycle. This approach revealed that the MEK-ERK cycle in PMA stimulated Jurkat-T cells could be operating in H or ST regimes depending on the strength of the stimulus [21]. A quasi-steady-state approximation Michealis-Menten model [22] employed in the hybrid approach [21] could not explain the observed transition of the regime effected by merely changing the stimulus strength. A question thus arises as to what could be the mechanism that may govern the observed operating regime transition under steady-state conditions. When Raf/MEK/ERK enzymatic cascades is modelled without an explicit feedback, signal flow is usually described as a one-way communication, that is, going from upstream to downstream of the cascade [23]. However, recently, a new type of signalling called retroactive signalling, which is caused by the presence of a downstream load, has been considered [24][25][26][27][28]. This phenomenon occurs due to the possibility that the futile cycles are coupled with another downstream cascade/substrate. Either or both forms of the protein involved in a futile cycle could be sequestered by another substrate which could be a part of another cascade or simply by a DNA to which one of the forms of the protein is sequestered [24-26, 29, 30]. In the case of Raf/MEK/ERK, the sequestration of the phosphorylated ERK could result in a retroactivity in the cascade. It has been shown experimentally that retroactivity indeed plays a role in the behaviour of MAPK cascades and other signalling pathways [24,[31][32][33][34][35]. Presence of retroactivity in enzymatic cascades has been suggested to predict a more realistic drug-response curve, that is, an input-output behaviour [30].
Inclusion of sequestration effects, which is known to affect the enzymatic futile cycle behaviour [24,35] may cause a shift in the operating regimes at deterministic level [30]. It is thus likely that incorporating the presence of retroactive signalling might predict the stimulus-strength dependent operating regime transition. In this study, we consider systematically characterising the effect of the presence of substrate or product retroactivity on the operating regimes of MEK/ERK enzymatic cascade. Specifically, we show that strength of the retroactive signalling can modulate the nature of the operating regimes and can permit operating regime transitions.

Mathematical model of a futile cycle with retroactive signalling
We consider an enzymatic futile cycle with retroactive signalling wherein an enzyme catalysed transition between inactive (M) and active (M p ) forms of the protein substrate occurs (Fig 1). We assume that both forms of the proteins M and M p may be sequestered reversibly, respectively by downstream targets S 1 and S 2 and thereby incorporate retroactivity in the cascade (Fig 1).
The biochemical reactions corresponding to the enzymatic cascade in Fig 1 are and those capturing the downstream sequestration steps are We assume quasi-steady state approximation (QSSA) for the two intermediate complexes in Eq (1) and (2), and for the two complexes formed by sequestration reactions (Eqs 3 and 4). Upon employing QSSA, the dynamics of dimensionless concentration of M p , � m = m p /m t where m t is the total protein substrate, dictated by the biochemical reactions in Eqs (1-4) is given by the mathematical kinetic model where, k f , and k r , respectively are the forward and reverse catalytic rate constants, e t and p t , respectively capture the total concentrations of kinase E and phosphatase P. K 1 and K 2 are the Michaelis-Menten (MM) constants for the forward and backward enzymatic reactions, respectively [22,36]. R p (e t , λ, � m) and R d (α = 0, � m), respectively capture the phosphorylation and dephosphorylation reaction rates. Assuming the equilibrium constants for binding of M and M p are equal, the retroactivity strengths for sequestration of M and M p , respectively are given by where s 1 and s 2 are the concentrations of species S 1 and S 2 , respectively and K d is the equilibrium constant corresponding to the sequestration reactions. A detailed derivation of Eq (5) from the full model capturing the dynamics of the biochemical reactions (Eqs 1-4) along with the definition of associated MM constants is in S1 Appendix. Note that the effect of retroactivity of either M or M p or both on the phosphorylation (first term in the right hand side or rhs) and dephosphorylation (second term in rhs) rates in Eq (5) is quantitatively accounted for by scaling the MM constants K 1 and K 2 with non-zero (positive) values of λ and α, respectively. Upon setting the left hand side or lhs to zero and solving analytically the resulting quadratic equation, we find the steady-state solution of Eq (5) as ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi b 2 À 4ðk f e t =k r p t Þð1 À k f e t =k r p t ÞðK 2 ð1 þ aÞ=m t Þ q 2ðk f e t =k r p t À 1Þ ; k f e t k r p t 6 ¼ 1 In S1 Text, we show that this steady-state solution (Eq 7) matches with that of the full model (Eq [AI.1-AI.4] and [AI.5]), for all range of values assigned to the parameters. Dose-response curve � m p ( � K 1 , � K 2 , e t ) of the futile cycle with (or without) retroactivity is essentially the locus of the relationship between � m p and e t , with all other parameters fixed [19,21]. Note that � m p can be drawn using Eq (7) for (a) without retroactivity by setting α = λ = 0, (b) with retroactivity only in M by setting α = 0, λ > 0, (c) with retroactivity only in M p by setting α>0, λ = 0 and (d) with retroactivity in both M and M p by setting α>0, λ>0 [24]. Since introduction of retroactivity tantamount to proportional scaling of the MM constants (Eq 5), for the sake of brevity, we define effective MM constants � K 1 = K 1 (1+λ) and � K 2 = K 2 (1+ α) which when λ or α set to zero will correspond to the case of absence of retroactivity in M or M p , respectively. Dose-response curve � m p can be classified into four distinct operating regimes, viz., H, ST, TH and U. Each of these regimes have a representative dose-response curve referred to as nominal profile. These four nominal profiles correspond to the four combinations of the saturated or unsaturated states of the two enzymatic reactions, viz., phosphorylation and dephosphorylation reactions of the futile cycle, as summarized in Table 1. An enzymatic reaction is considered saturated when most of the enzyme is bound to the substrate. The saturated state of the reaction occurs when the corresponding Michaelis-Menten constant is significantly smaller than the substrate concentration. � K n 1 , � K n 2 , where superscript n = H, ST, TH and U, used for arriving at the four nominal profiles of the futile cycle are in Table 1. As a ready reckoner, we present in  besides � K n 1 ; � K n 2 used for arriving at these curves are k r = k f = 0.01s -1 , p t = 200nM and m t = 1000nM [18,21]. Unless otherwise explicitly stated, these parameter values specified are employed for the rest of the study.
While hyperbolic response of the futile cycle is robust to fluctuations and can transmit signals in a broad range of amplitudes [37], signal transducing (ST) regime exhibiting a linear response is amenable for signalling involving graded stimuli. In the threshold-hyperbolic regime, the response of the futile cycle occurs only if the input is above the threshold, after which it increases hyperbolically [38]. Ultrasensitive (U) regime permits amplification of a small signal near the threshold which biological systems take advantage of [39]. As a reference, we employ the nominal profiles corresponding to the case wherein retroactivity is absent.
A dose-response curve is placed in one of the four regimes by contrasting the correspond- where superscript n = H, ST, TH, U indicates regime-specific nominal profile (Methods). This approach has been suggested by Gomez-Uribe et al. [18] and adopted in several recent studies [21,40].

Retroactivity impacts operating regimes
In order to study the effect of retroactivity on the dose-response curve, we adopt the same strategy prescribed by Gomez-Uribe et al. [18] to characterize the operating regimes in the presence of a downstream load on M or M p . We limit the scope of this study to the presence of retroactivity in either M or M p . Systematic characterization reported here, without loss of . The conditions employed for simulating these dose-response curves are in Table 1. In order to assess if retroactivity impacts the nature of the operating regime for a certain set of parameters, we consider a dose-response curve in the U regime in the absence of retroactivity (α = λ = 0), when ( � K 1 (0), � K 2 (0)) = (K 1, K 2 ) = (7,70).  Fig 3A). We further show that introduction of (a) retroactivity in M p can induce regime transition from H at �

Retroactivity strength dictates nature of regime transition
Retroactivity introduces a scaling for the Michaelis-Menten constants (Eq 5) and thereby affects the steady-state behaviour (Eq 7). As a result, in order to study the effect of retroactivity strength on the operating regimes, it is sufficient to understand how the parameter space of effective Michaelis-Menten constants � K 1 = K 1 (1+λ) and � K 2 = K 2 (1+α) is partitioned into different input-output behaviours. Note that replacing K 1 (1+λ) and K 2 (1+α) in Eq (7) respectively with � K 1 and � K 2 makes the retroactivity embedded steady-state solution form similar to that of an isolated enzymatic cascade. Thus, knowledge of the boundaries of the different operating regimes in the planes of � K 1 and � K 2 could be directly used to decipher the effect of  retroactivity on the dose-response curves exhibiting a certain input-output characteristic by varying λ or α. Next, we implemented an optimization problem to delineate the parameter space ( � K 1 , � K 2 ) corresponding to the four distinct operating regimes. For the ease of constructing the map, assuming α = λ = 0, for an operating regime, after specifying a � K 1 we identified � K 2 by increasing retroactivity strength α such that the candidate dose-response curve � m c p ð � K 1 ; � K 2 ; e t Þ satisfied the relative distance criterion ; e t Þ refers to a candidate.) This criterion is based on the metric suggested by Gomez-Uribe et al. [18] and recently used in Parundekar et al. [21]. In the metric introduced by Gomez-Uribe et al. [Gomez-Uribe et al.], the total substrate concentration is used as scaling. The predictions are therefore a function of the total substrate concentration itself. However, the basis for finding the distance from the nominal curves introduced by Parundekar et al. [21] constitutes scaling using the maximum regime-specific distance from its nominal profile. This metric offers advantages such as scaling being a self-learned parameter, relative distance estimation that is not biased by the system parameters. Next, we briefly describe the procedure adopted for estimating d c ( � K 1 , � K 2 , n). Every candidate dose-response curve will have four distances, each corresponding to a comparison with four regime-specific nominal profiles (Fig 2). Finding d c ( � K 1 , � K 2 , n) objectively for a dose-response curve requires estimation of max 8c k � m c p À � m n p k in Eq (8) a priori. However, the information about the regime to which a candidate � m c p belongs to is unavailable. In order to address this, we first created a randomly chosen parameter-profile database containing 140000 sets of ( � K 1 , � K 2 ) sampled using stratified random sampling (Methods) across five orders of magnitude range each tagged to its dose-response curve � m p . (Note that the maximum possible value that an element in � m p can take is 1 [21].) Next, we performed an optimization (Methods) for finding � K 2 that satisfies Eq (8) and its corresponding � m p . As an example, consider finding the boundaries of H regime by setting n = H in Eq (8). In the five-orders of magnitude range considered, finding ( � K 1 ; � K 2 ) whose corresponding � m p satisfied Eq (8) enabled identifying the boundary for the H regime in the planes of effective Michaelis-Menten constants (Fig 4, orange lines). Note that the dashed lines correspond to those ( � K 1 , � K 2 ) on the boundary sourced directly from the database. We repeated the entire procedure to find the boundaries corresponding to U (blue), ST (yellow), and TH (red) regimes (Fig 4). We note that upper boundary of the ST regime is an exception. While constructing the upper boundary for ST regime, we observed that the doseresponse curve is insensitive to � K 2 beyond a certain limit after which � K 2 has no effect on d c ( � K 1 , � K 2 , ST). Therefore, for representation purposes, we fixed the upper boundary for ST (yellow) at d c ð � K 1 ; � K 2 ; STÞ � 0:02 by accordingly modifying Eq 8. Note that as a direct consequence the dose-response curves well beyond the upper boundary of ST will belong to the signal-transducing regime. Metric adopted in Eq 8 by and large separates the regions where these four regimes exist. We note that the underlying model assumptions and the metric used by Gomez-Uribe et al. [18] are different as compared to those considered here. These differences could be attributed to the range for the operating regimes in Fig 4 not being same as those reported in [18].
Since changing retroactivity strength can independently modulate the Michaelis-Menten constants, manipulating �  (Fig 4, point D) when � K 2 is scaled to 60000. For a given source profile specified by a certain ( � K 1 , � K 2 ) with no retroactivity either in M or M p , while maintaining � K 1 or � K 2 constant, the minimum load λ min or α min , respectively required for inducing a regime transition is sensitive to the chosen � K 1 (0) or � K 2 (0) (S1 Text). This sensitivity analysis showed the minimum load needed for any regime transition to occur is 0.3. This minimum corresponds to transition of ST at ( � K 1 (λ = 0), � K 2 (α = 0)) = (8205,28160) to H regime due to retroactivity in M with λ min being 0.3. λ min = 0.3 translates to (ms 1 /(m tm p )) = λ min /(1+ λ min ) = 0.23 indicating that 23% of the unphosphorylated substrate sequestered by downstream target is needed for inducing this transition.

Saturation level of the two enzymatic reactions governs the retroactivity induced regime transition
Since the dose-response curve � m p explicitly depends on � K 1 (λ) and � K 2 (α) (Eq 7), understanding how load strength λ or α influences the input-output behaviour may offer useful insights into what causes retroactivity driven operating regime transition. In order to assess the regime-specific impact of retroactivity on the input-output behaviour, we systematically analyse the dose-response curves and the associated sensitivity with respect to retroactivity strengths λ and α.
The sensitivity of � m with respect to retroactivity strength λ and α, respectively are quantitatively captured by for a finite (non-zero) downstream load. Detailed expressions of these are in S2 Appendix. Eqs 9 and 10 show that the presence of retroactivity in M or M p introduces a constant scaling of � K 1 ð0Þ = K 1 or � K 2 (0) = K 2 , respectively to the sensitivity with respect to the corresponding Michaelis-Menten constant. In the sub-sections below, we present the sensitivity effects due to modulation of retroactivity corresponding to either M or M p for these five transitions and distil out the underlying causal mechanism. For the case of substrate or product retroactivity modulation, we first fix ( � K 1 ð0Þ = K 1 , � K 2 ð0Þ = K 2 ) in a certain regime with no retroactivity and then increasing λ or α, respectively and track the ensuing regime transition.

U to TH transition due to retroactivity in M.
The dose-response curves obtained by starting from U regime for ( � K 1 (λ = 0) = K 1 , � K 2 (α = 0) = K 2 ) = (9nM,9nM) with d c (9nM,9nM, U) = 0.013 and transitioning into TH regime by changing λ is shown in Fig 5A. Note that while α = 0, increasing λ leads to a proportional scaling of � K 1 (λ). Introduction of retroactivity causes changes to the extent of ultrasensitive nature of the dose-response curves. This extent of ultrasensitive nature in the presence of retroactivity can be quantified via the half-maximal response given by which uniquely specifies the dose-response curve's EC50, that is, e t at which � m = 0.5 [24]. Note that when λ = α = 0, Eq (11) reduces to the response defined in Goldbeter and Koshland [19]. As the dose-response curve transits from U to TH, the EC50 increases from 200 to 4400 for the range of λ considered (S1.3 Fig in S1 Text). Note that EC50 increases linearly with the retroactivity strength λ (Eq 11). Moreover, Fig 5A also reveals that an increase in load shifts the doseresponse curve by simultaneously enlarging the curve's base resulting in a threshold and also the curvature eventually leading to a TH input-output behavior. Next, we elucidate what causes the observed U to TH transition.
In Fig 5B, we show the modulation of sensitivity (Eq 9) by λ and e t . An increase in λ in dose-response curve � m p ( � K 1 (λ), � K 2 ð0Þ, e t ) leads to a decreased negative sensitivity. This shift in peak is correlated to the corresponding increase in the EC50 (S1.3 Fig in S1 Text), as has also been reported in Ventura et al. [24]. This behaviour is dictated by the steady-state levels of M p (Eq 7), which is a balance between the phosphorylation and dephosphorylation rate terms in the rhs of Eq 6 for a given λ and e t . Insights into the effect of λ on � m can be deciphered from the nature of relative variation of these two rates, which we discuss next.
In Fig 5C, we present the rate-balance plot consisting of the rate curves of the phosphorylation reaction R p (e t = 1000 nM, λ, � m) for different λ and of the dephosphorylation reaction R d (α = 0, � m). Note that R p and R d are the rates of the two enzymatic reactions defined in Eq (5). The nature of an enzymatic reaction being saturated, that is, all enzymes bound to its substrate, is specified by the range of � m for which the corresponding rate does not change significantly. Therefore a sufficient proportional increase in � K 1 (λ) due to λ can lead to R p exhibiting a linear dependence on � m in a certain range. The nature of the phosphorylation reaction is unsaturated in this range of � m. In the U regime, both R p (1000 nM, λ = 0, � m) (Fig 5C, black) and R d (0, � m) (Fig 5C, dashed) curves are predominantly saturated. For the chosen e t , at λ = 0, the intersection occurs in the region where R p is not saturated and R d is saturated, leading to � m � 1. Increasing λ forces the phosphorylation reaction (R p curve) to gradually become predominantly unsaturated (Fig 5C). The extent of this unsaturation introduced underlies the shift in the intersection point of the rate curves in the direction of decreasing � m. Thus, increasing λ causes significant decrease in the steady-state levels � m (Fig 5C). This decrease explains the gradual change in the steady-state levels at e t = 1000 nM in the different dose-response curves in Fig 5A. Moreover, this decrease results in a significant change in the sensitivity ( Fig  5B). At λ = 960, due to sufficient levels of unsaturation, the operating regime transits into the TH regime, which is characterized by the phosphorylation and dephosphorylation reactions, respectively being unsaturated and saturated. At the transition, the relative distance from the TH nominal profile � m TH p is d c ( � K 1 (960) = 8649, � K 2 = 9,TH) = 0.0833.

U to ST transition due to retroactivity in M p .
For ( � K 1 (λ = 0) = K 1 , � K 2 (α = 0) = K 2 ) = (7nM,70nM), the effect of dose-response curves on the retroactivity strength α is in Fig 5D. The dose-response curve when α = 0 (Fig 5D, black) with a d t (7nM,70nM, U) = 0.049 and EC50 of 178 nM, at α = 9.44 shifted to the ST regime with a d t ð � K 1 ð0Þ ¼ 7; � K 2 ð9:44Þ � 731; STÞ ¼ 0:0965, with the EC50 being 82 nM (S1.3 Fig in S1 Text). We next discuss what causes this regime transition. Fig 5E shows that an increase in the retroactivity strength α in dose-response curve � m p ( � K 1 (0), � K 2 (α), e t ) while maintaining λ = 0 causes reduction in the (positive) sensitivity. The ratebalance analysis at e t = 150nM shows that when α = 0, the intersection of the two rate-curves occurs in the region where the phosphorylation reaction is near saturation (Fig 5F). Note that the R d curve is predominantly saturated when α = 0. An increase in α, that is, � K 2 (α) shifts the nature of R d curve to predominantly unsaturated. This shift causes the intersection, that is, steady-state level, to increase from 0.2 at α = 0 to 0.99 at α = 18. Therefore, increasing the load leads to an increase in the steady-state level depending on the extent of the unsaturation evidenced by the dephosphorylation reaction. This shift in the steady-state level forces the doseresponse curve to move into the ST operating regime.

H to ST transition due to retroactivity Mp.
For ( � K 1 (λ = 0) = K 1 , � K 2 (α = 0) = K 2 ) = (3000nM,3000nM), the effect of dose-response curves on the retroactivity strength α is in Fig  5G. At α = 0, the dose-response curve belonged to the H regime with a d c (3000nM,3000nM, H) = 0.014. Upon increasing α to 5, dose-response curve transitions to ST operating regime with d c ð � K 1 ð0Þ ¼ 3000; � K 2 ð5Þ � 18000; STÞ ¼ 0:08. EC50 for the dose-response curves changes from 200nM to~38nM (S1.3 Fig in S1 Text). Increase in the load causes a decrease in the sensitivity. The rate-balance plot for e t = 150nM shows that both phosphorylation and dephosphorylation reactions are predominantly unsaturated in the H regime for α = 0. Upon increasing the load, while the dephosphorylation reaction continues to remain unsaturated, the rate curve shows a slowed-down response to increase in � m, that is, reduction of the slope of the R d curve. This reduction causes a shift in the intersection of the rate-curves to a larger substrate concentration. For e.g., for e t = 150nM, the steady-state level at α = 0 and 20 are 0.42 and 0.94, respectively. Thus, the dose-response curve transitioning from the H to ST is essentially caused by this reduction in the slope of the R d curve with increase in the retroactivity strength in M p .

ST to H and TH to H transitions due to retroactivity M and Mp.
In Fig 6, we show the dose-response curves capturing the regime transition from ST to H and TH to H driven by load in M and M p , respectively. For the case of ST to H transition (Fig 6A), at ( � K 1 (λ = 0) = K 1 , � K 2 (α = 0) = K 2 ) = (9 nM,9000nM), the dose-response curve has an EC50 of~11 with a d c (9,9000,ST) = 0.003. With an increase in the retroactivity strength λ in M, the EC50 increases and at λ = 400, the dose-response curves achieved corresponds to the H regime with a d c ( � K 1 (400) = 3609, � K 2 = 9000,TH) = 0.092 and EC50 of 86.5 (S1.3 Fig in S1 Text). Sensitivity analysis and rate-balance plot at e t = 12nM show that while dephosphorylation reaction is unsaturated, an increase in λ the phosphorylation reaction transitions from predominantly saturated to unsaturated state and thereby, driving the regime transition (S1.4 Fig in S1 Text).

Discussion and conclusion
Input-output behaviour of an activated enzymatic futile cycle has been studied extensively due to its ability to orchestrate cell fate in direct and indirect context-dependent manners [2,14,21]. Michealis-Menten constants (MM) dictated saturated/unsaturated state of the two enzymatic reactions facilitates placing steady-state dose-response curves of a futile cycle into Signal transducing (ST), Hyperbolic (H), Threshold hyperbolic (TH) and Ultrasensitive (U) operating regimes (Fig 2) [18]. The unphosphorylated (M) and phosphorylated (M p ) forms of the protein substrate involved in the futile cycle can be sequestered by respective downstream

PLOS ONE
targets. The sequestration dictated load or retroactivity on the upstream protein levels introduces a two-way signal flow permitting modulation of the steady-state behaviour [24,26]. In this study, we systematically show that the presence of retroactivity in M or M p can shift the input-output behaviour from one operating regime to another by modulating the level of saturation or unsaturation of the enzymatic reactions. In particular, we demonstrate five possible transitions: (a) U to TH and ST to H caused by retroactivity in M and (b) U to ST, TH to H, and H to ST by that in M p .
Using a quasi-steady state approximated model of the futile cycle with retroactivity, we systematically identified the MM constants' range that permit four distinct operating regimes in the presence of retroactive signalling (Fig 4). Surprisingly, the minimum retroactivity strength needed for inducing any transition is 0.3 which translates to 23% of the substrate bound to its target. For this minimum retroactivity strength of 0.3, dose-response curve at ( � K 1 (λ = 0), � K 2 (α = 0)) = (8205,28160) belonging to ST regime transitions into H regime. Several downstream targets that could sequester proteins in MAPK cascades have been reported [35]. Thus, while analysing such a behaviour in a futile cycle using experimental data, ignoring the hidden retroactive signalling effect, however small, could lead to an incorrect prediction of the underlying operating regime.
While in this study we only considered increasing the retroactivity strength to trigger a regime transition, in principle, if downstream sequestrations were already present, its strength can be decreased too. Decreasing the retroactivity strength can predict five other transitions that are essentially the reverse of those analysed in this study. Further, simultaneous increase (decrease) of the retroactivity strengths in M and M p can lead to a transition from U to H (H to U). Thus, the operating regime boundaries reported in Fig 4 permits prediction of all 12 possible regime transitions. We further note that the U and H regimes have a slight overlap in the ( � K 1 , � K 2 ) space. Our study indicates that recent experimental observations that a stimulusstrength dependent shift in the operating regimes is possible in a single MAPK cascade if the stimulus concentration change could cause retroactivity induced regime transition [21].
Using sensitivity and rate-balance analysis, we demonstrated that modulation of the saturation or unsaturation levels of the two enzymatic reactions by changing the retroactivity strength is the fundamental reason for the operating regime transition. In particular, we show that increase in retroactivity (a) in M leads to increasing unsaturation in the phosphorylation reaction and (b) in M p makes dephosphorylation reaction more unsaturated (Fig 5, S1.4 and S1.5 Figs in S1 Text). This is due to the fact that the steady-state level of M p , the active form, is sensitive to changes in the retroactivity strength. While increasing the strength of retroactivity in M causes a decrease in the (negative) sensitivity of the steady-state level, that in M p leads to marked reduction in the (positive) sensitivity. This sensitivity to retroactive signalling can be capitalized upon to modulate the nature of response of the futile cycle. Synthetic biology tools are becoming available for tweaking the binding sites of targets to which the protein substrate, active/inactive forms may bind and thereby enabling control of the extent of sequestration [41,42]. The nature of sensitivity effect that retroactive signalling bestows on the steady-state levels demonstrated in this study can be of immense value for precise engineering of a cell to control and modulate the input-output behaviour.

Regime identification
The regime that a dose-response curve � m p belongs to is identified by contrasting it with the four nominal profiles � m n p ð � K n 1 ; � K n 2 ; e t Þ where superscript n = H, ST, TH, U. The dose-response curve � m p is ascribed to a certain regime H, ST, TH, or U if the relative distance between � m p and the corresponding nominal profile is within 10%.

Stratified random sampling
The two stratification cut-off points were chosen. While choosing these cut-off points, (a) 60000 samples were chosen in the (0<K 1 <1600, 0<K 2 <1600) and (b) 10000 samples each in the range (0<K 1 <50, 0<K 2 <10000) and (0<K 1 <10000, 0<K 2 <50) were chosen [21]. In both these cases uniform distribution was used for sampling. Samples corresponding to either one or both reactions being saturated was at least 10% more than that for the case where both reactions were unsaturated.

Optimization for finding operating regime boundaries
The boundary for a specific regime was obtained by seeking � K 2 that satisfied the objective function in Eq (8) solved using nonlinear optimization function "fmincon" implemented in Matlab 1 [43]. A tolerance of 1e -6 was set as convergence criteria to the optimization problem. Optimizer convergence was sluggish in the presence of steep gradients and in these cases, the ( � K 1 , � K 2 ) samples from the database that satisfied Eq (8) was used.
Supporting information S1 Text. S1.1-Steady-state solution of the full model. S1.2-Dose-response curve transition from Signal-Transducing to Hyperbolic regime induced by substrate retroactivity. S1.3-Minimum retroactivity strength required to transition from one regime to another. S1.4-Effect of retroactivity strength on EC50 during different regime transitions. S1.5-Sensitivity and rate balance analysis to investigate retroactivity induced ST to H and TH to H regime transitions.