Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Retroactivity induced operating regime transition in an enzymatic futile cycle

  • Akshay Parundekar,

    Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing – review & editing

    Affiliation Department of Chemical Engineering, Indian Institute of Technology Bombay, Powai, Mumbai, India

    ⨯
  • Ganesh A. Viswanathan

    Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing

    ganeshav@iitb.ac.in

    Affiliation Department of Chemical Engineering, Indian Institute of Technology Bombay, Powai, Mumbai, India

    ⨯

Abstract

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.

Citation: Parundekar A, Viswanathan GA (2021) Retroactivity induced operating regime transition in an enzymatic futile cycle. PLoS ONE 16(4): e0250830. https://doi.org/10.1371/journal.pone.0250830

Editor: Christopher Rao, University of Illinois at Urbana-Champaign, UNITED STATES

Received: November 11, 2020; Accepted: April 14, 2021; Published: April 30, 2021

Copyright: © 2021 Parundekar, Viswanathan. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the manuscript and its Supporting Information files.

Funding: This study was supported by Science and Engineering Research Board, Department of Science and Technology, Government of India (MTR/2020/000589 and CRG/2020/002672) for funding this study. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

1. 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–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–12]. Aberrant functioning of this cascade has been implicated in many diseases such as cancer [13, 14]. Detailed understanding of the sustained and transient activation patterns of 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–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–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–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.

2. 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 (Mp) forms of the protein substrate occurs (Fig 1). We assume that both forms of the proteins M and Mp may be sequestered reversibly, respectively by downstream targets S1 and S2 and thereby incorporate retroactivity in the cascade (Fig 1).

thumbnail
Download:
Fig 1. Enzymatic futile cycle with retroactivity.

M and Mp are the inactive and active forms of the protein substrate. Kinase E and phosphatase P, respectively are the enzymes for the phosphorylation and dephosphorylation biochemical reactions. While S1 and S2 are the downstream targets, respectively of M and Mp, MS1 and MpS2 are the corresponding sequestered complexes.

https://doi.org/10.1371/journal.pone.0250830.g001

The biochemical reactions corresponding to the enzymatic cascade in Fig 1 are [1] [2] and those capturing the downstream sequestration steps are [3] [4]

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 Mp, = mp/mt where mt is the total protein substrate, dictated by the biochemical reactions in Eqs (1–4) is given by the mathematical kinetic model [5] where, kf, and kr, respectively are the forward and reverse catalytic rate constants, et and pt, respectively capture the total concentrations of kinase E and phosphatase P. K1 and K2 are the Michaelis-Menten (MM) constants for the forward and backward enzymatic reactions, respectively [22, 36]. Rp(et, λ, ) and Rd(α = 0, ), respectively capture the phosphorylation and dephosphorylation reaction rates. Assuming the equilibrium constants for binding of M and Mp are equal, the retroactivity strengths for sequestration of M and Mp, respectively are given by [6] where s1 and s2 are the concentrations of species S1 and S2, respectively and Kd 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 Mp 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 K1 and K2 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 [7] where, . 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 (, , et) of the futile cycle with (or without) retroactivity is essentially the locus of the relationship between and et, with all other parameters fixed [19, 21]. Note that 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 Mp by setting α>0, λ = 0 and (d) with retroactivity in both M and Mp 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 = K1 (1+λ) and = K2 (1+ α) which when λ or α set to zero will correspond to the case of absence of retroactivity in M or Mp, respectively. Dose-response curve 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. , , 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 Fig 2 the nominal dose-response curves with n = H, ST, TH, U for the four regimes. Parameters besides used for arriving at these curves are kr = kf = 0.01s-1, pt = 200nM and mt = 1000nM [18, 21]. Unless otherwise explicitly stated, these parameter values specified are employed for the rest of the study.

thumbnail
Download:
Fig 2.

Schematic showing the steady-state dose response curve corresponding to the nominal profiles of (A) Hyperbolic (), (B) Signal transducing (), (C) Threshold hyperbolic (), (D) Ultrasensitive (). The conditions employed for simulating these dose-response curves are in Table 1.

https://doi.org/10.1371/journal.pone.0250830.g002

thumbnail
Download:
Table 1. The nature of the state of the two biochemical reactions corresponding to four operating regimes.

https://doi.org/10.1371/journal.pone.0250830.t001

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 corresponding (, , et) with , 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].

3. Results

3.1 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 Mp. We limit the scope of this study to the presence of retroactivity in either M or Mp. Systematic characterization reported here, without loss of generalization, can be used for the case where retroactivity may be present in both M and Mp, simultaneously.

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 ( (0), (0)) = (K1, K2) = (7,70). Fig 3A shows this dose-response curve (solid yellow) contrasted against the nominal profile for U regime (dashed blue), included from Fig 2D for ease of comparison, used for identifying the regime to which it belongs to. Introduction of retroactivity in Mp with a strength of α = 27 (and λ = 0) resulting in ( (0), (27) = (7,1960)) causes shifting of the dose-response curve (solid purple curve in Fig 3A) to the left. for case of α = 27 belongs to the ST regime indicating the possibility of retroactivity induced transition of operating regimes. (For the sake of comparison, we present (dashed red curve) from Fig 2B in Fig 3A). We further show that introduction of (a) retroactivity in Mp can induce regime transition from H at (0), (0) = (389,751) to ST regime at ((0), (1.5) = (389,1878) (Fig 3B) and (b) retroactivity in M can induce operating regime transition from ST to H (S1 Text). Given that the presence of a downstream load can cause a regime shift, we ask a question as to what are the other possible transitions in the presence of retroactivity. The primary goal of this study is to systematically understand the effect of retroactivity in M or Mp on the operating regimes.

thumbnail
Download:
Fig 3.

Retroactivity in Mp inducing operating regime transition from (A) U at (λ = 0) = 7, (α = 0) = 70 to ST at ((0), (27)) = (7,1960) and (B) H at ((λ = 0) = 389, (α = 0) = 751) to ST at ( (0), (15)) = (389,12016). Inset: Zoom in of the dose-response curves. For ease of comparison, the nominal profiles and from Fig 2D and 2B, respectively are included in (A). Similarly, and from Fig 2A and 2B, respectively are included in (B). Parameters ((0), (0)) used for simulating the nominal profiles are in Table 1.

https://doi.org/10.1371/journal.pone.0250830.g003

3.2 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 = K1(1+λ) and = K2(1+α) is partitioned into different input-output behaviours. Note that replacing K1(1+λ) and K2(1+α) in Eq (7) respectively with and 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 and 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 (, ) corresponding to the four distinct operating regimes. For the ease of constructing the map, assuming α = λ = 0, for an operating regime, after specifying a we identified by increasing retroactivity strength α such that the candidate dose-response curve satisfied the relative distance criterion [8] for all n = H, ST, TH, and U. (Note that superscript c in 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 dc(, , n).

Every candidate dose-response curve will have four distances, each corresponding to a comparison with four regime-specific nominal profiles (Fig 2). Finding dc(, , n) objectively for a dose-response curve requires estimation of in Eq (8) a priori. However, the information about the regime to which a candidate belongs to is unavailable. In order to address this, we first created a randomly chosen parameter-profile database containing 140000 sets of (, ) sampled using stratified random sampling (Methods) across five orders of magnitude range each tagged to its dose-response curve . (Note that the maximum possible value that an element in can take is 1 [21].) Next, we performed an optimization (Methods) for finding that satisfies Eq (8) and its corresponding . 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 () whose corresponding 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 (, ) 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 dose-response curve is insensitive to beyond a certain limit after which has no effect on dc(, , ST). Therefore, for representation purposes, we fixed the upper boundary for ST (yellow) at 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].

thumbnail
Download:
Fig 4. Boundaries of the four operating regimes hyperbolic (H, blue), signal transducing (ST, green), threshold-hyperbolic (TH, red), and ultrasensitive (U, yellow) in the planes of and for λ = α = 0.

All boundaries of each of the regimes except the upper boundary of ST satisfy the relative distance criterion in Eq (8). For the case of upper boundary of ST the rhs of Eq (11) was set to 0.02. (, ) on the dotted lines extending the solid line boundaries were sourced directly from the database. While dose-response curves corresponding to parameter sets at green dots A and B were used as example of transition from U to TH regime, those curves at C and D of transition from H to ST.

https://doi.org/10.1371/journal.pone.0250830.g004

Since changing retroactivity strength can independently modulate the Michaelis-Menten constants, manipulating or or both could cause a shift in the characteristic input-output behaviour. Specifically, by increasing the strength of the load in M causing proportional change in while keeping constant, a dose-response curve in U or ST, respectively can shift to TH or H. For example, (10,10) in U regime (Fig 4, point A) would shift to TH (Fig 4, point B) upon increasing to 10000. Similarly, while maintaining constant, an increase in the retroactivity strength in Mp leading to proportional change in could lead to four other possible regime transitions, viz., U to ST, TH to H or ST, and H to ST. (6000,6000) in H regime (Fig 4, point C) transitions into ST regime (Fig 4, point D) when is scaled to 60000. For a given source profile specified by a certain (, ) with no retroactivity either in M or Mp, while maintaining or constant, the minimum load λmin or αmin, respectively required for inducing a regime transition is sensitive to the chosen (0) or (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 ((λ = 0), (α = 0)) = (8205,28160) to H regime due to retroactivity in M with λmin being 0.3. λmin = 0.3 translates to (ms1/(mt—mp)) = λmin/(1+ λmin) = 0.23 indicating that 23% of the unphosphorylated substrate sequestered by downstream target is needed for inducing this transition.

3.3 Saturation level of the two enzymatic reactions governs the retroactivity induced regime transition

Since the dose-response curve explicitly depends on (λ) and (α) (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 with respect to retroactivity strength λ and α, respectively are quantitatively captured by [9] and [10] 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 Mp introduces a constant scaling of = K1 or (0) = K2, 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 Mp for these five transitions and distil out the underlying causal mechanism. For the case of substrate or product retroactivity modulation, we first fix ( = K1, = K2) in a certain regime with no retroactivity and then increasing λ or α, respectively and track the ensuing regime transition.

3.3.1 U to TH transition due to retroactivity in M.

The dose-response curves obtained by starting from U regime for ((λ = 0) = K1, (α = 0) = K2) = (9nM,9nM) with dc(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 (λ). 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 [11] which uniquely specifies the dose-response curve’s EC50, that is, et at which = 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 dose-response 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.

thumbnail
Download:
Fig 5.

Effect of retroactivity strength on the operating regimes and the associated sensitivity for (i) U to TH, (ii) U to ST, and (iii) H to ST transitions. While panel (i) corresponds to effects due to load on M quantified by λ, panel (ii) and (iii) captures those due to load on Mp quantified by α. Dependence of dose-response curves on the load corresponding to (i), (ii) and (iii) are in (A), (D) and (G), respectively. Sensitivity of steady-state level for different retroactivity strengths for (i), (ii), and (iii) are in (B), (E) and (H), respectively. Sensitivity curves in (B) was estimated using Eq (9), Eq (10) was employed for those in (E) and (H). While rate-balance plot showing the effect of retroactivity strength on modulation of steady-state levels corresponding to (i) at et = 1000nM is in (C), that for (ii) and (iii) at et = 150nM are in (F) and (I), respectively. Colorbar in each of the panels display the retroactivity strengths. Dotted line in (A), (B) and (C) in panel (i) corresponds to the dose-response, sensitivity and Rd(λ) curves, respectively at the transition where λ = 960. Dotted line in panels (ii) and (iii) captured these curves at the corresponding transition where α = 9.44 and α = 5, respectively.

https://doi.org/10.1371/journal.pone.0250830.g005

In Fig 5B, we show the modulation of sensitivity (Eq 9) by λ and et. An increase in λ in dose-response curve ((λ), , et) 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 Mp (Eq 7), which is a balance between the phosphorylation and dephosphorylation rate terms in the rhs of Eq 6 for a given λ and et. Insights into the effect of λ on 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 Rp(et = 1000 nM, λ, ) for different λ and of the dephosphorylation reaction Rd(α = 0, ). Note that Rp and Rd 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 for which the corresponding rate does not change significantly. Therefore a sufficient proportional increase in (λ) due to λ can lead to Rp exhibiting a linear dependence on in a certain range. The nature of the phosphorylation reaction is unsaturated in this range of . In the U regime, both Rp(1000 nM, λ = 0, ) (Fig 5C, black) and Rd (0, ) (Fig 5C, dashed) curves are predominantly saturated. For the chosen et, at λ = 0, the intersection occurs in the region where Rp is not saturated and Rd is saturated, leading to . Increasing λ forces the phosphorylation reaction (Rp 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 . Thus, increasing λ causes significant decrease in the steady-state levels (Fig 5C). This decrease explains the gradual change in the steady-state levels at et = 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 is dc((960) = 8649, = 9,TH) = 0.0833.

3.3.2 U to ST transition due to retroactivity in Mp.

For ((λ = 0) = K1, (α = 0) = K2) = (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 dt(7nM,70nM, U) = 0.049 and EC50 of 178 nM, at α = 9.44 shifted to the ST regime with a , 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 ( (0), (α), et) while maintaining λ = 0 causes reduction in the (positive) sensitivity. The rate-balance analysis at et = 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 Rd curve is predominantly saturated when α = 0. An increase in α, that is, (α) shifts the nature of Rd 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 dose-response curve to move into the ST operating regime.

3.3.3 H to ST transition due to retroactivity Mp.

For ((λ = 0) = K1, (α = 0) = K2) = (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 dc(3000nM,3000nM, H) = 0.014. Upon increasing α to 5, dose-response curve transitions to ST operating regime with . 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 et = 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 , that is, reduction of the slope of the Rd curve. This reduction causes a shift in the intersection of the rate-curves to a larger substrate concentration. For e.g., for et = 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 Rd curve with increase in the retroactivity strength in Mp.

3.3.4 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 Mp, respectively. For the case of ST to H transition (Fig 6A), at ((λ = 0) = K1, (α = 0) = K2) = (9 nM,9000nM), the dose-response curve has an EC50 of ~11 with a dc(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 dc((400) = 3609, = 9000,TH) = 0.092 and EC50 of 86.5 (S1.3 Fig in S1 Text). Sensitivity analysis and rate-balance plot at et = 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).

thumbnail
Download:
Fig 6.

Dose-response curves capturing the retroactivity driven transition of operating regimes from (A) ST to H and (B) TH to H. Colorbar displays the retroactivity strength corresponding to the dose-response curves. Dotted lines in (A) and (B), respectively correspond to the retroactivity strength λ = 400 and α = 578 at which the regime transition occurs.

https://doi.org/10.1371/journal.pone.0250830.g006

The dose-response curve at ((λ = 0) = K1, (α = 0) = K2) = (9000 nM,9 nM) belonging to the TH regime with dc(9000,9,TH) = 0.05 transitions into H regime with increase in load α from 0 to 578 (Fig 6B). The dc((0) = 9000, (578) = 5211,H) = 0.07. Sensitivity and rate-balance plots at et = 3000 nM suggests that regime transition is caused by the phosphorylation reaction being unsaturated and the dephosphorylation reaction shifting from being predominantly saturated at α = 0 to primarily unsaturated with increasing α (S1.5 Fig in S1 Text).

4. 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 (Mp) forms of the protein substrate involved in the futile cycle can be sequestered by respective downstream 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 Mp 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 Mp.

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 ((λ = 0), (α = 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 Mp 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 (, ) space. Our study indicates that recent experimental observations that a stimulus-strength 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 Mp 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 Mp, 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 Mp 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.

5. Methods

5.1 Regime identification

The regime that a dose-response curve belongs to is identified by contrasting it with the four nominal profiles where superscript n = H, ST, TH, U. The dose-response curve is ascribed to a certain regime H, ST, TH, or U if the relative distance between and the corresponding nominal profile is within 10%.

5.2 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<K1<1600, 0<K2<1600) and (b) 10000 samples each in the range (0<K1<50, 0<K2<10000) and (0<K1<10000, 0<K2<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.

5.3 Optimization for finding operating regime boundaries

The boundary for a specific regime was obtained by seeking that satisfied the objective function in Eq (8) solved using nonlinear optimization function “fmincon” implemented in Matlab® [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 (, ) 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.

https://doi.org/10.1371/journal.pone.0250830.s001

(PDF)

S1 Graphical abstract.

https://doi.org/10.1371/journal.pone.0250830.s002

(TIFF)

S1 Appendix. Derivation of the mathematical kinetic model.

https://doi.org/10.1371/journal.pone.0250830.s003

(DOCX)

S2 Appendix. Sensitivity of steady-state level to retroactivity strength.

https://doi.org/10.1371/journal.pone.0250830.s004

(DOCX)

Acknowledgments

We acknowledge the DBT-Pan IIT Center for Bioenergy, IIT Bombay for an access to the high-performance computing facility.

References

  1. 1. Widmann C, Gibson S, Jarpe MB, Johnson GL. Mitogen-activated protein kinase: conservation of a three-kinase module from yeast to human. Physiol Rev. 1999;79(1):143–80. pmid:9922370
  2. 2. Chang L, Karin M. Mammalian MAP kinase signalling cascades. Nature. 2001;410(6824):37–40. pmid:11242034
  3. 3. Saurin AT. Kinase and Phosphatase Cross-Talk at the Kinetochore. Frontiers in Cell and Developmental Biology. Vol 6, 2018. p. 62. pmid:29971233
  4. 4. Kleiman LB, Maiwald T, Conzelmann H, Lauffenburger DA, Sorger PK. Rapid phospho-turnover by receptor tyrosine kinases impacts downstream signaling and drug binding. Mol Cell. 2011 Sep;43(5):723–37. pmid:21884975
  5. 5. Blüthgen N, Bruggeman FJ, Legewie S, Herzel H, Westerhoff H V, Kholodenko BN. Effects of sequestration on signal transduction cascades. FEBS J. 2006;273(5):895–906. pmid:16478465
  6. 6. Ortega F, Acerenza L, Westerhoff H V, Mas F, Cascante M. Product dependence and bifunctionality compromise the ultrasensitivity of signal transduction cascades. Proc Natl Acad Sci. 2002;99(3):1170–5. pmid:11830657
  7. 7. Thattai M, Van Oudenaarden A. Intrinsic noise in gene regulatory networks. Proc Natl Acad Sci. 2001;98(15):8614–9. pmid:11438714
  8. 8. Dhananjaneyulu V, Nanda Sagar VP, Kumar G, Viswanathan GA. Noise Propagation in Two-Step Series MAPK Cascade. PLoS One. 2012;7(5):35958. pmid:22563473
  9. 9. Baraskar AA, Deb A, Viswanathan GA. Noise Propagation in Series Enzymatic Cascades. IFAC Proc Vol. 2013;46(31):89–94.
  10. 10. Ondrey FG, Dong G, Sunwoo J, Chen Z, Wolf JS, Crowl-Bancroft C V, et al. Constitutive activation of transcription factors NF-(kappa)B, AP-1, and NF-IL6 in human head and neck squamous cell carcinoma cell lines that express pro-inflammatory and pro-angiogenic cytokines. Mol Carcinog. 1999 Oct;26(2):119–29. pmid:10506755
  11. 11. Aggarwal S, Takada Y, Singh S, Myers JN, Aggarwal BB. Inhibition of growth and survival of human head and neck squamous cell carcinoma cells by curcumin via modulation of nuclear factor‐κB signaling. Int J cancer. 2004;111(5):679–92. pmid:15252836
  12. 12. Yang J, Li G, Zhang K. Pro-survival effects by NF-κB, Akt and ERK(1/2) and anti-apoptosis actions by Six1 disrupt apoptotic functions of TRAIL-Dr4/5 pathway in ovarian cancer. Biomed Pharmacother. 2016 Dec;84:1078–87. pmid:27780136
  13. 13. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144(5):646–74. pmid:21376230
  14. 14. Dhanasekaran DN, Johnson GL. MAPKs: function, regulation, role in cancer and therapeutic targeting. Oncogene 2007 May 14;26(22):3097–9. pmid:17496908
  15. 15. Santos SDM, Verveer PJ, Bastiaens PIH. Growth factor-induced MAPK network topology shapes Erk response determining PC-12 cell fate. Nat Cell Biol 2007 Mar 18;9(3):324–30. pmid:17310240
  16. 16. Frank SA. Input-output relations in biological systems: measurement, information and the Hill equation. Biol Direct. 2013;8(1):31. pmid:24308849
  17. 17. Sepulchre J-A, Merajver SD, Ventura AC. Retroactive signaling in short signaling pathways. PLoS One. 2012;7(7):e40806. pmid:22848403
  18. 18. Gomez-Uribe C, Verghese GC, Mirny L a. Operating regimes of signaling cycles: statics, dynamics, and noise filtering. PLoS Comput Biol 2007 Dec;3(12):e246.
  19. 19. Goldbeter A, Koshland DE Jr. An amplified sensitivity arising from covalent modification in biological systems. Proc Natl Acad Sci U S A 1981 Nov;78(11):6840–4. pmid:6947258
  20. 20. Fey D, Halasz M, Dreidax D, Kennedy SP, Hastings JF, Rauch N, et al. Signaling pathway models as biomarkers: Patient-specific simulations of JNK activity predict the survival of neuroblastoma patients. Sci Signal. 2015;8(408):ra130–ra130. pmid:26696630
  21. 21. Parundekar A, Kalantre G, Khadpekar A, Viswanathan GA. Operating regimes in a single enzymatic cascade at ensemble-level. PLoS One. 2019;14(8).
  22. 22. Michaelis L, Menten ML. aDie Kinetik der Invertinwirkung, o Biochem. z. 1913;49:333–69 (Johnson KA, Goody RS. The original Michaelis constant: translation of the 1913 Michaelis–Menten paper. Biochemistry. 2011;50(39):8264–9)
  23. 23. Kholodenko BN, Hancock JF, Kolch W. Signalling ballet in space and time. Nat Rev Mol Cell Biol 2010 Jun;11(6):414–26. pmid:20495582
  24. 24. Ventura AC, Jiang P, Van Wassenhove L, Del Vecchio D, Merajver SD, Ninfa AJ. Signaling properties of a covalent modification cycle are altered by a downstream target. Proc Natl Acad Sci U S A. 2010/05/17. 2010 Jun 1;107(22):10032–7. pmid:20479260
  25. 25. Ventura AC, Sepulchre J-A, Merajver SD. A hidden feedback in signaling cascades is revealed. PLoS Comput Biol. 2008 Mar;4(3):e1000041. pmid:18369431
  26. 26. Del Vecchio D, Ninfa AJ, Sontag ED. Modular cell biology: retroactivity and insulation. Mol Syst Biol 2008/02/12. 2008;4:161. pmid:18277378
  27. 27. Catozzi S, Di-Bella JP, Ventura AC, Sepulchre J-A. Signaling cascades transmit information downstream and upstream but unlikely simultaneously. BMC Syst Biol. 2016;10(1):84. pmid:27561377
  28. 28. Grunberg TW, Del Vecchio D. Modular analysis and design of biological circuits. Curr Opin Biotechnol. 2020;63:41–7. pmid:31877425
  29. 29. Ossareh HR, Ventura AC, Merajver SD, Del Vecchio D. Long signaling cascades tend to attenuate retroactivity. Biophys J 2011 Apr 6;100(7):1617–26. pmid:21463574
  30. 30. Ventura AC, Jackson TL, Merajver SD. On the role of cell signaling models in cancer research. Cancer Res. 2009;69(2):400–2. pmid:19147549
  31. 31. Jiang P, Ventura AC, Sontag ED, Merajver SD, Ninfa AJ, Del Vecchio D. Load-induced modulation of signal transduction networks. Sci Signal. 2011;4(194):ra67–ra67. pmid:21990429
  32. 32. Jayanthi S, Nilgiriwala KS, Del Vecchio D. Retroactivity controls the temporal dynamics of gene transcription. ACS Synth Biol. 2013;2(8):431–41. pmid:23654274
  33. 33. Brewster RC, Weinert FM, Garcia HG, Song D, Rydenfelt M, Phillips R. The transcription factor titration effect dictates level of gene expression. Cell. 2014;156(6):1312–23. pmid:24612990
  34. 34. Kim Y, Coppey M, Grossman R, Ajuria L, Jiménez G, Paroush Z, et al. MAPK substrate competition integrates patterning signals in the Drosophila embryo. Curr Biol. 2010 Mar 9;20(5):446–51. pmid:20171100
  35. 35. Kim Y, Paroush Z, Nairz K, Hafen E, Jiménez G, Shvartsman SY. Substrate-dependent control of MAPK phosphorylation in vivo. Mol Syst Biol 2011 Feb 1;7:467. pmid:21283143
  36. 36. Puskas J.E., Kantor J., Shrikande G., 2017. Reaction engineering with enzymes: A relatively uncharted territory. AIChE J. 63,266–272.
  37. 37. Levine J, Kueh HY, Mirny LA. Intrinsic fluctuations, robustness and tunability in signalling. Biophys J. 2007. 92, 4473–4481.
  38. 38. Birtwistle M, Rauch J, Kiyatkin A, Aksamitine E, Doberzynski M, Hoek J, et al. Emergence of bimodal cell population responses from the interplay between analog single-cell signaling and protein expression noise. BMC Syst Biol. 2012. 6, 109. pmid:22920937
  39. 39. Jr Ferrell. JE, Machleder EN. The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes. Science. 1998. 280, 895–898.
  40. 40. Straube R. Operating regimes of covalent modification cycles at high enzyme concentrations. J Theor Biol. 2017;431:39–48. pmid:28782551
  41. 41. Morsut L, Roybal KT, Xiong X, Gordley RM, Coyle SM, Thomson M, et al. Engineering Customized Cell Sensing and Response Behaviors Using Synthetic Notch Receptors. Cell 2016 Feb 11;164(4):780–91. pmid:26830878
  42. 42. Zhou X.X., Fan L.Z., Li P., Shen K., Lin M.Z., Optical control of cell signaling by single-chain photoswitchable kinases. Science. 2017. 355, 836–842. pmid:28232577
  43. 43. Mathworks, Inc. 2017. Weblink: http://in.mathworks.com/help/optim/ug/fmincon.html