Frequency-preference response in covalent modification cycles under substrate sequestration conditions

Abstract

Covalent modification cycles (CMCs) are basic units of signaling systems and their properties are well understood. However, their behavior has been mostly characterized in situations where the substrate is in excess over the modifying enzymes. Experimental data on protein abundance suggest that the enzymes and their target proteins are present in comparable concentrations, leading to substrate sequestration by the enzymes. In this enzyme-in-excess regime, CMCs have been shown to exhibit signal termination, the ability of the product to return to a stationary value lower than its peak in response to constant stimulation, while this stimulation is still active, with possible implications for the ability of systems to adapt to environmental inputs. We characterize the conditions leading to signal termination in CMCs in the enzyme-in-excess regime. We also demonstrate that this behavior leads to a preferred frequency response (band-pass filters) when the cycle is subjected to periodic stimulation, whereas the literature reports that CMCs investigated so far behave as low-pass filters. We characterize the relationship between signal termination and the preferred frequency response to periodic inputs and we explore the dynamic mechanism underlying these phenomena. Finally, we describe how the behavior of CMCs is reflected in similar types of responses in the cascades of which they are part. Evidence of protein abundance in vivo shows that enzymes and substrates are present in comparable concentrations, thus suggesting that signal termination and frequency-preference response to periodic inputs are also important dynamic features of cell signaling systems, which have been overlooked.

Introduction

Biological systems must respond to internal and external variations such as the depletion of nutrients, the fluctuations in hormone levels, and the arrival of sensory signals. In response to stimuli, the pathway-controlling enzymes change their activities. Two basic phenomena play a significant role in this processing: allosteric changes in protein conformation and covalent modification of proteins¹.

Covalent modification cycles (CMCs) are one of the major intracellular signaling mechanisms, both in prokaryotic and eukaryotic organisms². In such cycles, a signaling protein is modified by the addition of a chemical group. This modification may result in either activation or inactivation, depending on the particular signaling pathway involved, followed by a reverse process, thus closing the cycle. For phosphorylation–dephosphorylation cycles, two opposing enzymes are involved: a kinase and a phosphatase. In the absence of external stimulation, the cycle is in a steady-state where the activation and inactivation reactions are dynamically balanced. External stimuli that produce a change in the enzymatic activity shift the activation state of the target protein, creating a departure from a steady state, which can propagate through a signaling cascade. While individual CMCs are simply elements of a large signaling network, understanding their response to inputs is an essential first step in characterizing the response of more-elaborated signaling networks to external stimuli.

There is a large body of literature on CMCs and cascades of CMCs in the substrate-in-excess regime using mathematical modeling tools^1,3,4,5,6,7. This regime implies that the substrate abundance is in large excess over the modifying enzymes’ abundances. Importantly, it was predicted that these systems can be highly sensitive to changes in stimuli if their catalyzing enzymes are saturated with their target protein substrates¹ (note that enzyme saturation is a stronger requirement than substrate in excess). This mechanism was termed zero-order ultrasensitivity and has received enormous attention throughout the years.

However, the substrate-in-excess condition cannot be guaranteed in vivo, because endogenous enzyme concentrations are much higher than those used in a typical in vitro assay^8,9. In cascades, experimental data on protein abundance suggest that both enzymes and their target proteins are present in comparable concentrations. This is the case of the mitogen-activated protein kinase (MAPK) cascade, which was studied using a combination of theoretical and experimental tools¹⁰. Estimation of the parameters associated with each level of the cascade demonstrated that the amount of protein in both the second and third levels are similar. In general, in a cascade, the protein at one level operates as the enzyme in the next one, suggesting that at least the two levels involved are not within the substrate-in-excess regime. Enzyme-in-excess scenarios (the opposite of substrate-in-excess) were studied recently supported by in vivo levels of kinases and phosphatases frequently exceeding the levels of their corresponding substrates in budding yeast¹¹.

Theoretical work focusing on the departure from the substrate-in-excess plus enzyme saturation conditions is scarce^12,13. In particular, it is not clear how CMCs in an enzyme-in-excess regime and the cascades of which these CMCs are partly responding to external inputs. It was shown that sequestration of the substrate results in a reduction in ultrasensitivity, which changes the dynamics of a CMC and may account for signal termination and a sign-sensitive delay¹⁴, in which the rise in the signal is delayed but the dropping signal is transduced immediately. Such a delay element provides cells with units that neglect short fluctuations in signals but transduce long signals¹⁵. In another study, the importance of sequestration-based feedback in signaling cascades was thoroughly analyzed, and a positive feedback mechanism that emerges from sequestration effects was shown to bring about bistability in the cascade¹⁶. Negative feedback has been shown to emerge from sequestration effects, as was theoretically predicted⁷ and then experimentally validated^17,18. Furthermore, the input–output curves were classified in terms of the saturation state of the activating and inactivating enzymes, including the ultrasensitive regime³. Only two of the mentioned studies have addressed the response properties of CMCs to fluctuating external inputs^3,17.

In this paper, we investigate the response properties of CMCs in the enzyme-in-excess regime to both step-constant stimulations (inputs abruptly increasing from zero to a value that remains constant in time) and periodic stimulation resulting from sequestration of the substrate protein. The complexity of the transient response patterns of a dynamical system to constant stimulation ranges from a monotonic increase (relatively simple), to overshoot (intermediate) to damped oscillations. They reflect the different effective time scales present in the system, which are uncovered by the input and result from the interplay of the system’s time constants, and may have different functional consequences. Overshoot responses reflect a property of the underlying system referred to as adaptation or signal termination; i.e., their ability to return to a steady-state value lower than the peak response while the stimulation is still active. Requirements for biochemical adaptation were extensively studied^19,20. Additional scenarios leading to the same type of signal termination responses include receptor internalization and receptor desensitization²¹ and the response to protein variation²². To our knowledge, there are no reports of signal termination in the substrate-in-excess regime, whereas, in the enzyme-in-excess regime, signal termination was predicted¹⁴, but not investigated in detail. The questions arise of whether signal termination is present only in this enzyme-in-excess regime and under what conditions this occurs. The study of relatively simple systems shows that signal termination may arise from negative feedback or incoherent feedback loop mechanisms¹⁹. Therefore, understanding the mechanisms underlying signal termination will shed light on the broader problem of characterizing the role of negative regulators in cell signaling²³. Previous work on cascades in the substrate-in-excess regime has shown the presence of damped oscillations in CMC cascades⁷ in response to constant inputs, raising the question of whether these patterns persist in the enzyme-in-excess regime.

The transient dynamics associated with adaptation reflect an effective time scale of the system, which was shown to have implications for their preferred frequency response to periodic inputs^24,25,26,27. In spite of this, the frequency response of CMCs has received much less attention than adaptive behavior. Some studies found that CMCs behave as tunable low-pass filters, filtering out high-frequency fluctuations or noise in signals and environmental cues^3,17,28,29. A frequency-preference response was reported in CMCs only under a dose-conservation scheme (dose-conservation implies that either by amplitude or duration compensation, the total dose is kept constant when the frequency increases)³⁰. However, all these studies are away from the enzyme-in-excess regime and thus, they preclude signal termination in the unforced CMC.

In this paper, we use mathematical modeling and detailed computational simulations to characterize the conditions under which CMCs in the enzyme-in-excess regime exhibit adaptation to constant inputs (signal termination) and preferred frequency responses to periodic inputs, their properties, how these two phenomena are related, and their consequences for the dynamic behavior of the cascades these CMCs are part of. Specifically, we focus on the regimes where enzymes are in similar or higher amounts than the substrates, resulting in a large fraction of the substrate protein sequestrated by the enzymes. Cascades are not simply chains of CMCs, where one CMC feeds the subsequent one, but because of the backward connections between CMCs, cascades are relatively complex networks⁷ and therefore their dynamics cannot be simply inferred from the dynamics of the CMC components.

Our modeling approach involves the so-called mechanistic models fully describing the interactions between enzymes, substrates, and products overall the time scales present in the process. The mathematical descriptions of the functioning of a CMC usually reduce the system’s dimensionality by means of a quasi-steady-state approximation (QSSA)³¹, leading to a one-dimensional system. However, one-dimensional systems do not exhibit signal termination and preferred frequency responses to periodic inputs in realistic conditions (i.e., unless they are imposed to the system), and therefore these approximations are not applicable to the CMCs we study. Subsequent approximations applicable to the substrate-in-excess regime (Goldbeter & Koshland, 1981) or the so-called total QSSA (tQSSA), which is applicable also when the concentrations of substrate and enzyme are comparable or the enzyme-in-excess regimes (Tzafriri, 2003), also lead to one-dimensional systems and therefore they are also not applicable to our study either. The main reason for these failures is that the fast time scale, which is neglected in the above-mentioned dimensionality process, plays a significant role in producing the two phenomena that are the object of our study. To validate these ideas, we compare our results using the “full” models with the results using approximated models.

Understanding how a response is triggered is as important as deciphering why it terminates and how it can be reactivated. Different time scales emerge from these processes, the response time, the duration of the signal, the recovery time. If the stimulating signal has an associated natural time scale as well, as is the case for time-varying signals, it is expected that interesting behavior could emerge from the interaction of all the involved time scales. Given that in vivo data suggest that sequestration of the substrate is at play in CMCs and in signaling pathways where those CMCs are involved, and given that those pathways are subjected to time-varying signals, understanding their frequency response in sequestration conditions is not only interesting but also relevant.

To summarize, the main focus of the paper is on frequency-preference response in CMCs under substrate sequestration conditions. The paper deals with two system-level behaviors that are not well characterized in CMCs: signal termination and frequency-preference response. It also deals with the enzyme-in-excess regime, a regime that, while physiological, was not properly studied in the literature. The paper studies the connection between having a transient response under constant stimulation, referred to as signal termination, and exhibiting a frequency-preference response under periodic stimulation. Finally, the article also evaluates how the modeling approximations that reduce systems’ dimensionality fail in describing the temporal properties of the responses. These last two features (connection transient response—frequency-preference response, and failure of some simplifying approximations) can be extrapolated to other systems and are not exclusive to CMCs.

The organization of the paper is as follows. We first study CMCs with a mechanistic mathematical model and under substrate sequestration conditions, characterizing signal termination in the parameter space. We then apply periodic stimulation to CMCs exhibiting signal termination, finding conditions for frequency-preference emergence. We then study cascades of CMCs. Finally, we evaluate the performance of the approximations usually employed to describe CMCs, in detecting both signal termination and frequency preferences. We conclude by discussing the potential implications of the results in the article.

Results

Signal termination in CMCs under sequestration conditions

To our knowledge, the only report on signal termination in CMCs focused on the interplay of two enzymes, kinase (activating) and phosphatase (inactivating) in 3:1 concentrations relative to the substrate, thus leading to sequestration conditions¹⁴. The CMC exhibited a transient (overshoot) response under constant stimulation, therefore, terminating the prolonged response signal while the stimulation is still active. This required a fast kinase with low affinity and a slow phosphatase with high affinity. The fast kinase phosphorylates the available target, but the phosphorylated target is subsequently sequestered by the low-activity high-affinity phosphatase. In the steady-state, most of the target substrate is sequestered by phosphatase. In this section, we investigate whether the three conditions that were found to ensure signal termination in CMCs (enzymes in excess over the substrate, kinase faster than phosphatase, and phosphatase with higher affinity than kinase) are necessary conditions and their relative importance. For the remaining of this article, we refer to kinase and phosphatase as the activating and inactivating enzymes, respectively. Our results are general and apply to CMCs of the form presented in Fig. 1a.

Fig. 1: Signal termination in a CMC: the importance of the relative concentrations between substrate and enzymes.

a Scheme representing a CMC cycle. S is the unmodified substrate, P the product (modified substrate), E₁ and E₂ are the modifying enzymes (kinase and phosphatase), C₁ is the substrate–kinase complex, C₂ is the product-phosphatase complex. Three kinetic parameters are associated with each enzymatic reaction, a_i, d_i, and k_i: the association, dissociation, and catalytic rates, respectively (more details in the Methods section). b Phosphorylated protein (P, solid green line), substrate–kinase complex (C₁, dotted pink line), and product-phosphatase complex (C₂, dashed blue line) are plotted versus time. Parameter values: a₁ = 3.5; d₁ = 1; k₁ = 50; a₂ = 0.3; d₂ = 0.25; k₂ = 0.25; E_1T = 100; E_2T = 100; S_T = 35. a_i, i = 1,2, are in units of 1/(time × concentration), d_i and k_i are in units of 1/time, the total enzymes (E_1T and E_2T) are in units of concentration (see Methods section). P_max and P_ss, the maximum and the steady-state levels of P, are indicated over the plot. The characteristic value 0.63 × P_max is indicated with an asterisk over the time course. c Graphs of Q (=P_max/P_ss) versus E_1T/S_T and E_2T/S_T, with Q in color scale. E_1T/S_T, E_2T/S_T, and Q are in logarithmic scale, E_1T/S_T and E_2T/S_T vary between 0.1 and 10, resulting in Q values in the range between 1 and ~46. The solid blue line over the plots corresponds to Q = 1.6 and separates signal termination from monotonic behavior and non-monotonic behavior not satisfying the criteria. The solid black line over the plot corresponds to P_ss = 0.2. The cases over the line (P_ss > 0.2) are termed strong signal termination. Dashed black lines correspond to P_ss = 0.1 or 0.3, indicating different boundaries for the strong signal termination region. On the upper left triangle, Q < 1.6 (no signal termination satisfying the criteria). The kinetic parameter values are as in b. d Representative time courses in different regions in c. Parameter set 1: E_1T/S_T = 0.13, E_2T/S_T = 0.13, parameter set 2: E_1T/S_T = 1, E_2T/S_T = 0_.13, parameter set 3: E_1T/S_T = 1, E_2T/S_T = 0_.79, parameter set 4: E_1T/S_T =1, E_2T/S_T = 1_.26, parameter set 5: E_1T/S_T = 0.13, E_2T/S_T = 7.9, parameter set 6: E_1T/S_T = 7.9, E_2T/S_T = 7.9.

In Fig. 1, we illustrate the phenomenon of signal termination in CMCs (Fig. 1a) for a representative set of kinetic parameter values and analyze the importance of the relative concentrations of the substrate and enzymes (results with an alternative parameter set are included in the Supplementary Information, Supplementary Figure 7). The total amount of substrate (S_T), kinase (E_1T), and phosphatase (E_2T) are conserved. S_T is set to zero for times previous to t = 0 s and follows a step-like profile for t > 0 s. Studies with step-like profiles starting at non-zero levels are included in the Supplementary Information (Supplementary Figure 9). Similar step-like stimulations but in different parameters are considered in the Supplemental Information (Supplementary Figures 1 and 2) and a discussion on why considering S_T as a signal is included in the Discussion. Additional details and the description of the mechanistic mathematical model used here are presented in Methods. Figure 1b shows the time course of P, C₁, and C₂ under similar conditions as in¹⁴. C₁ and P increase fast (C₁ faster than P) and C₂ increases on a much slower time scale. The signal (P) terminates (reaches its maximum P_max and then relaxes to its steady-state value P_ss), whereas the stimulation is still active, because of the sequestration of P by C₂.

In what follows, it is reasonable to focus on signals that exhibit significant termination according to some criteria, which we define here. First, we normalize the time courses of each variable by their conserved total amount, so that each variable is within 0 and 1, as in Fig. 1b. We neglect those outputs for which P_max < 0.1. Second, we select those outputs for which P_ss < 0.63 × P_max. Finally, we require P_ss < 0.2. This ensures a strong enough signal with a large enough decay to steady-state, which in turn ensures a clear signal peak, and a return to values of P close to the pre-stimulation values. We refer to the signals satisfying these three conditions as strong signal termination. For future use, we define the peak-to-steady state amplitude Q = P_max/P_ss. The second condition implies Q > 1.5873 (~1.6). Although the specific numbers chosen (0.1, 0.63, and 0.2) are somehow arbitrary, similar results to the ones presented in this paper are obtained for other choices (Supplementary Section 5).

In Fig. 1c, we evaluate how signal termination depends on the relative concentrations of kinase/substrate (E_1T/S_T) and phosphate/substrate (E_2T/S_T). To this end, we scanned E_1T and E_2T within some range (lower, equal, and higher than S_T), whereas keeping all other (kinetic) parameters fixed (same values as in Fig. 1b). The blue curve corresponds to Q = 1.6 (log(Q) = 0.47)) and separates regions for which Q > 1.6 (above, satisfying the decay condition) and Q < 1.6 (below). Points above this line correspond to signal termination. A subset of these points (in between the two solid black curves) exhibits strong signal termination (signals satisfy the three conditions discussed in the previous paragraph). On the left upper triangle, there is no signal termination (neither strong nor mild) and the outputs are such that P_max < 0.1. Examples of the signal (P) behavior in each region are presented in Fig. 1d.

Figure 1c demonstrates that the occurrence of both signal termination and strong signal termination is primarily controlled by E_2T/S_T. This means that the phosphatase has to be in (roughly) similar or higher concentrations than the substrate for the two behaviors to occur, while there are weaker restrictions on the relative concentrations of the kinase and substrate. Supplementary Figure 7 in the Supplementary Information shows a similar analysis for a different choice of kinetic parameters.

We next characterize the emergence of signal termination and its dependence on the kinetic parameters of the CMC model. We will focus on the two kinetic conditions mentioned above (kinase faster than phosphatase, phosphatase with higher affinity than the kinase), which are satisfied for the parameter values used in Fig. 1b (the kinase association and dissociation rates are one order of magnitude higher than those for the phosphatase and two orders higher in the case of the catalytic rate).

The velocity V of an enzymatic reaction is the rate at which the product is formed. For the two reactions in the CMC cycle, these are given by³¹:

$$V = \frac{{\rm{d}}{P}}{\rm{dt}} = V_{{\rm{max}},1}\frac{S}{{K_{m,1} + S}},\;V = \frac{{\rm{d}}{S}}{\rm{dt}} = V_{{\rm{max}},2}\frac{P}{{K_{m,2} + P}}$$

(1)

The associated parameters V_max and K_m for each reaction are, respectively, the maximum reaction velocity, defined as the catalytic rate of the enzyme multiplied by its total amount, and the Michaelis constant K_m = (d+k)/a that indicates the amount of substrate, leading to a velocity V_max/2 (we omit the indices of a, d, and k for clarity). Equation 1 indicates that at small substrate concentrations the rate of the reaction is linear, with an effective rate constant V_max/K_m, while at large substrate concentrations the reaction saturates at its maximum rate V_max. For a fixed amount of enzyme E_T, the analysis of the role of the velocity of each reaction in the emergence of signal termination reduces to the analysis of k/K_m, and k.

The affinity Aff = 1/K_m of an enzyme for the substrate measures the concentration of substrate that must be present to saturate the enzyme. A high (low) value of Aff indicates that a small (large) concentration of substrate is needed to saturate the enzyme.

Simple calculations using the parameter values used in Fig. 1b confirm that the kinase is faster than the phosphatase but has a lower affinity than the phosphatase. In fact, K_m1 = 14.6 and K_m2 = 1.7 (Aff₁ ~ 0.07 and Aff₂ ~ 0.59) and the two indicators of the reaction’s velocity, k_i/K_m,i and k_i, result in 3.4 and 50 for the kinase, 0.14 and 0.25 for the phosphatase.

In order to evaluate the impact of relative velocities and relative affinities of kinase and phosphatase in producing signal termination, we performed the following study where the corresponding model parameters for the opposing enzymes (a₁ and a₂, d₁ and d₂, k₁ and k₂) varied along lines in parameter space with slopes α_a, α_d, and α_k, respectively. More specifically, (a₁, d₁, k₁) = (α_a x a₂, α_d x d₂, α_k x k₂) in the range 0.1–10. When the two opposing enzymes have the same kinetics (α_a = 1, α_d = 1, α_k = 1), these kinetics were selected to be those of the phosphatase in Fig. 1b. This phosphatase-based study (kinase relative to the phosphatase) was complemented with an analogous kinase-based study (phosphatase relative to the kinase) using a parameter β with the same characteristics of the parameter α above, but when the two opposing enzymes have the same kinetics, these kinetics were selected to be those of the kinase in Fig. 1b.

For the phosphatase-based study, starting from that parameter set, we consider four cases: a common control parameter for the three kinetic parameters (case 1) and a control parameter for any one of the kinetic parameters, whereas the other two remain fixed and equal (cases 2, 3, 4). More specifically,

• 1. α_a = α_d = α_k = α, meaning (a₁, d₁, k₁) = α(a₂, d₂, k₂),

• 2. α_d = α_k = 1 and α_a = α, meaning a₁ = α a₂ and(d₁, k₁) = (d₂, k₂),

• 3. α_a = α_k = 1 and α_d = α, meaning d₁ = α d₂ and (a₁, k₁) = (a₂, k₂),

• 4. α_a = α_d = 1 and α_k = α, meaning k₁ = α k₂ and (a₁, d₁) = (a₂, d₂).

The description of the four cases considered for the kinase-based study is analogous.

For all these cases, we derive analytically the effect of the control parameters (α and β) over the relative effective first-order rate constants V₁/V₂ and the relative affinities Aff₂/Aff₁, accordingly, where V₁ = k₁/K_m,1 and V₂ = k₂/K_m,2 (see Supplementary Section 2). A simultaneous increase in V₁/V₂ and Aff₂/Aff₁ with increasing values of α favors signal termination, so it is expected that Q increases with α. The opposite effects on V₁/V₂ and Aff₂/Aff₁ act against signal termination, so it is expected that Q decreases with α. Any other situation produces competing effects, so the behavior of Q is not straightforward, this behavior indicates which effect has stronger control in signal termination. The resulting Q also depends on the reference parameter set.

Our results are presented in Fig. 2. The top panels correspond to the phosphatase-based study and the bottom panels correspond to the kinase-based study. The shadowed region corresponds to Q < 1.6 (no signal termination as discussed above). The separating line coincides with V₁=V₂ and Aff₁=Aff₂. In all cases, we describe changes as either α or β increase.

Fig. 2: Signal termination in a CMC: kinase faster than phosphatase, phosphatase with higher affinity.

a Velocity of the kinase (E₁) relative to that of the phosphatase (E₂) and affinity of the phosphatase relative to that of the kinase (phosphatase parameters values as in Fig. 1b). In each plot, V₁/V₂ and Aff₂/Aff₁ are indicated in the left axis while Q (=P_max/P_ss) is indicated in the right axis. V₁/V₂ is plotted in green, Aff₂/Aff₁ in purple, and Q in black. The shadowed region over the plot indicates values of Q < 1.6 (no signal termination). b Velocity and affinity of the phosphatase relative to those of the kinase (kinase parameters as in Fig. 1b). Other details are as in a. Parameter values: a₁ = 3.5; d₁ = 1; k₁ = 50; a₂ = 0.3; d₂ = 0.25; k₂ = 0.25; E_1T = 100; E_2T = 100; S_T = 35.

We first analyze the phosphatase-based cases (top panels). In Case 1 (Fig. 2a, Case 1), V₁/V₂ increases, whereas Aff₂/Aff₁ remains constant, resulting in Q increasing with α for α higher than 1. Changes along the a₁–a₂ line (Case 2, Fig. 2a) result in an increase in V₁/V₂ but a decrease of Aff₂/Aff₁, leading to competing effects: the kinase increases its velocity but also its affinity with α. The resulting Q is 1 for all values of α, indicating that in this scenario the responses are monotonically increasing. Changes along the d₁–d₂ line (Case 3, Fig. 2a) result in a decrease in V₁/V₂ and an increase in Aff₂/Aff₁, leading to competing effects: the kinase decreases its velocity and its affinity with α. As in Case 2, the resulting Q is 1 for all values of α. Finally, changes along with the k₁–k₂ line (Case 4, Fig. 2a) result in an increase in both V₁/V₂ and Aff₂/Aff₁. Both effects favor signal termination, resulting in an increase in Q versus α.

A similar analysis for the kinase-based cases shows that in Cases 1 and 2 (Fig. 2b) V₁/V₂ decreases, whereas Aff₂/Aff₁ is constant in Case 1 and increasing in Case 2. Case 1 results in signal termination for β < 0.5, whereas Case 2 does not present signal termination because the kinase is faster than the phosphatase for β < 1, whereas the phosphatase has a higher affinity for β > 1. Case 3 (Fig. 2b) exhibits constant values of V₁/V₂ and Aff₂/Aff₁, resulting in values of Q that do not exceed the threshold for signal termination. Case 4 (Fig. 2b) exhibits almost constant V₁/V₂ too, but in this case, Aff₂/Aff₁ is a decreasing function of β, resulting in signal termination for β < 0.6.

Summarizing the results in Fig. 2, for the two studies, the only two cases resulting in signal termination are Cases 1 and 4 (Fig. 2a, b). Case 4 satisfies the two requirements, V₁/V₂ > 1 and Aff₂/Aff₁ > 1, for α > 1 or for β < 1, so it was expected to achieve signal termination. In Case 1, in contrast, kinase and phosphate have the same affinity, but still results in signal termination when α or β are such that V₁/V₂ > 1. Cases with competing effects in the relations V₁/V₂ and Aff₂/Aff₁ resulted in no signal termination. All the calculations related to Fig. 2 are included in Supplementary Section 2.

We now extend the analysis discussed above, which is a local analysis around a selected parameter set (that of Fig. 1b) with the goal of extracting more general conclusions. To this end, we performed a random parameter space exploration in the following ranges: E_1T, E_2T, and S_T in 10–100, the association and dissociation rates (a₁, d₁, a₂, d₂) in 0.1–10, and the catalytic rates (k₁, k₂) in 1–50. These parameter ranges were determined after a numerical study of the system and based on the literature, as detailed in the Methods section. Modified ranges are presented in Supplementary Figure 8 and do not affect the conclusions of this section. For each simulation, we randomly selected the parameter values within these ranges following the description in Methods, using a Latin Hypercube Sampling. We classified the outputs according to whether they produce signal termination or not. Previous work, mentioned above, provided the following conditions for signal termination¹⁴: E_1T/S_T > 1, E_2T/S_T > 1, V₁/V₂ > 1, Aff₂/Aff₁ > 1. Therefore, we also classified the output according to whether these conditions are satisfied or not.

Our results are presented in Fig. 3. In Fig. 3a, b, the yellow dots note all the points in parameter space for which simulations were performed (20,000). The cyan dots note the points in parameter space for which Bluthgen’s conditions are satisfied (696, 3,48%). By construction, the “Bluthgen dots” are located in the first quadrant where two of the four conditions (E_1T/S_T > 1, E_2T/S_T > 1, V₁/V₂ > 1, Aff₂/Aff₁ > 1) are evidently satisfied. If the other two are satisfied as well, these dots in the first quadrant are colored cyan. The purple dots note the points in parameter space for which the output shows signal termination (according to the criteria discussed above) (561, 2.80% from which 442, 2.21%, show strong signal termination). Finally, black dots indicate the overlap between cyan and purple dots, i.e., parameter sets satisfying Bluthgen’s conditions and signal termination (254, 1.27%). The existence of purple not overlapping cyan dots indicates that Bluthgen’s conditions can be relaxed and still obtain signal termination. These conditions cannot be fully violated without losing the signal termination phenomena, as indicated by the fact that no signal termination points lie on the third quadrant of Fig. 3a, and only a few in Fig. 3b. Instead, the relaxation can occur either by extending the conditions to the second quadrant where Aff₂/Aff₁ < 1 or E_1T/S_T < 1, or to the fourth quadrant where V₁/V₂ < 1 or E_2T/S_T < 1. Only a few points exhibit signal termination with E_2T/S_T < 1 showing that it is possible, but not robust. Figure 3b also shows that signal termination is more sensitive to condition E_2T/S_T > 1 than to condition E_1T/S_T > 1: a significantly higher number of signal termination cases lie in the second quadrant than in the fourth one. There are only a few outlier cases not satisfying both conditions (only 4 out of 20,000 within this study).

Fig. 3: Parameter space exploration.

a V₁/V₂–Aff₂/Aff₁ parameter space, b E_2T/S_T–E_1T/S_T parameter space. Both panels include the output of numerical simulations for 20,000 parameter sets in the selected ranges of variation. Each dot represents a single simulation output. Dots in purple correspond to outputs with signal termination, dots in cyan correspond to the parameter sets that satisfy Bluthgen’s conditions (see text), dots in black correspond to the parameter sets that satisfy Bluthgen’s conditions and signal termination. All the other simulation results are in yellow. Outputs lower than 0.1 were excluded from the analysis. Since E_1T and E_2T are not kept constant (as was the case in Fig. 2), V_i = k_iE_1T/K_mi with i = 1,2. c V₁/V₂–Aff₂/Aff₁ parameter space, but plotting only outputs with signal termination and magnifying the signal termination region. Each output is color and symbol coded according to the value of E_2T/S_T and E_1T/S_T, gray dots if E_2T/S_T and E_1T/S_T > 1, cyan stars if only E_2T/S_T > 1, yellow triangles if only E_1T/S_T > 1, green crosses if E_2T/S_T and E_1T/S_T < 1. d E_2T/S_T−E_1T/S_T parameter space, but plotting only outputs with signal termination and magnifying the signal termination region. Each output is color and symbol coded according to the value of V₁/V₂ and Aff₂/Aff₁, gray dots if V₁/V₂ and Aff_2/Aff₁ > 1, cyan stars if only V₁/V₂ > 1, yellow triangles if only Aff₂/Aff₁, green crosses if V₁/V₂ and Aff₂/Aff₁ < 1. Horizontal and vertical black lines over the plots indicate V₁ = V₂, Aff₂ = Aff₁, E_1T = S_T, E_2T = S_T, respectively. The characteristics of each distribution are captured in box plots giving the median (color line as central value), the 95% confidence interval of the median, the first and third quartiles (box), the 5th, and 95th percentiles (end of whiskers). Extra symbols in line with the box plots are outliers.

Figure 3c presents the results of the same simulations, plotting only those cases showing signal termination with a magnification of the quadrants of interest. Figure 3c has a color and symbol code for cases with E_1T/S_T and E_2T/S_T > 1, only E_2T/S_T > 1, only E_1T/S_T > 1, none of them. Figure 3d has a color and symbol code for cases with V₁/V₂ and Aff₂/Aff₁ > 1, only Aff₂/Aff₁ > 1, only V₁/V₂ > 1, none of them. Median values for each group and in each coordinate are indicated with a boxplot on top and on the right side of each panel (details in the Methods section).

A comparative analysis of signal termination and strong signal termination outputs is included in Supplementary Figure 3. From that analysis, we conclude that low values of V₁/V₂ and high values of Aff₂/Aff₁ and of E_2T/S_T lead to the strong signal termination regime. With E_1T/S_T it is not possible to distinguish a region that clearly promotes strong signal termination.

Frequency-preference response of CMCs to periodic inputs under sequestration conditions

In this section, we study the response of CMCs to time-varying stimuli, in conditions such that they exhibit signal termination under constant stimulation. The presence of adaptation in autonomous linear (and some nonlinear) dynamical systems has been associated with their ability to exhibit resonance, a peak in their amplitude response to oscillatory inputs at a preferred (resonant) frequency. Similar band-pass filtering behavior was reported in other biochemical systems and owing to different mechanisms^32,33,34. In the case of CMCs, signal termination and frequency-preference responses result from the complex interaction of effective time scales where the slow time constant of the negative feedback has a prominent role. The question arises whether and under what conditions CMCs that exhibit signal termination also exhibit preferred responses to oscillatory inputs. Although we hypothesized the occurrence of resonant-like responses in CMCs, this is not obvious since the effective time scales leading to signal termination depends on a combination of model parameters (rates) that govern the dynamics of autonomous CMCs.

The impact of the input frequency variation on the response amplitude to oscillatory periodic inputs is typically evaluated by measuring the gain, defined as the output amplitude normalized by the input amplitude as a function of the input frequency. This definition is an extension of the standard impedance for linear circuits, applicable to nonlinear systems under a number of assumptions including (i) the number of input and output cycles coincide, and (ii) the steady-state output amplitude is uniform across cycles when the input oscillations have this property. For simplicity, this type of study is usually conducted by keeping the input amplitude constant for all input frequencies. However, as we discuss below, the presence of conservation laws conditions the ability to use this type of inputs without violating the non-negativity of the substrate (S) concentration, and adapted (modified) versions of the input must be used. We note that we are using the word “adaptation” to refer to signal termination in the context of the response to constant inputs and to the modifications needed for the periodic inputs as discussed here. The two are unrelated concepts.

More specifically, we measure the product (P) response to periodic fluctuations in the total substrate (S_T). We define the gain as the quotient of the amplitudes of P and S_T (see details below). In this section, we aim to apply a periodic control over S_T in a frequency-independent manner by using either sinusoidal and square waves (Fig. 4; expected S_T). We refer to this as the expected variation since, as we discuss here, it cannot always be achieved. As indicated in the conservation condition in the scheme of Fig. 1a, the substrate is not always free (S), but it is mainly sequestered by the phosphatase (to form the complex C₂) and also, but in a lower amount, forming the product (P), or forming the complex with the kinase (C₁). So, in most scenarios, it is not possible to directly control S_T to follow the expected variation. The closest possible way to do this involves controlling free S, either by adding it or washing it out. These changes have an impact on S_T and how it is distributed among the different species. However, when S reaches the zero level and the expected variation in S_T indicates that S_T has to decrease even further, the task cannot be achieved instantaneously. Instead, one must cease to add S and allow the reactions to progress until a new balance is reached leading to the generation of free S, which can be then washed out allowing a further decrease in S_T. However, this decrease in S_T due to the recently released S being washed out has the time scale given by the reactions, generating a waveform different from that for the expected variation. Furthermore, the minimum value reached by S_T cannot be arbitrarily selected, but it is determined by the process just described. This situation is illustrated in Fig. 4, for both the expected sinusoidal and square wave variation in S_T and their adapted versions, which we use in our simulations.

Fig. 4: Periodic control of the total amount of substrate.

a–c Sinusoidal variation in S_T. The expected variation in S_T (purple solid) is given by S_T = 35(1+0.5sin(ωt))/1.5. The accomplished variation in S_T (green dashed) is adapted as explained in the text. The scale on the right of the figure corresponds to S (free substrate, solid yellow). The horizontal axis starts at 50 s in order to include stationary variations only. The arrow (a) indicates the time at which complete depletion of S occurs (omitted in the other panels). a Parameter values: a₁ = 3.5; d₁ = 1; k₁ = 30; a₂ = 0.3; d₂ = 0.25; k₂ = 2.5; E_1T = 100; E_2T = 100; S_T = 35. The frequency of stimulation is f = 0.02 Hz, corresponding to T = 50 s. b Parameter values are as in a with the exception of a₁ and a₂ that are multiplied by 0.1, resulting in a₁ = 0.35 and a₂ = 0.03. In this way, the system is overall slower than in a. For a slower CMC, the differences between expected and adapted are greater. The frequency of stimulation is f = 0.02 Hz. c Parameter values as in a but with a frequency of stimulation of f = 0.04 Hz, corresponding to T = 25 s. The stimulation is faster than in a, resulting in greater differences between expected and adapted S_T. d–f Train of square pulses in S_T. CMC parameter values as in a–c.

Al alternative approach to adapt the periodic stimulation is evaluated in Supplementary Figure 6 together with an analytical calculation that clarifies the ideas behind the need of adapting the forcing S_T (from expected to adapt) due to the sequestration of some components appearing in the conservation law. This calculation is done in a scenario that is simpler than a CMC, where analytical calculations are not possible (Supplementary Section 4).

We define the optimal frequency response of the CMC to variations in S_T as the ability of the amplitude or the gain to peak at a non-zero finite frequency (Fig. 5b, d). We measure the amplitudes of S_T and P as the difference between the maximum and minimum of these quantities once the output has reached the stationary, periodic variation regime (i.e., disregarding the transients). The amplitude and gain profiles (curves of the output amplitude and gain as a function of the input frequency) are presented in Fig. 5a, b (left: amplitude profiles, right: gain profiles). For both types of input, sinusoidal (Fig. 5a) and square waves (Fig. 5b), the CMC exhibits optimal frequency responses in the gain (Fig. 5b, c, right) respectively, for parameter values for which the underlying autonomous system exhibits signal termination. Figure 5b, d (left) illustrates that the output amplitude does not necessarily capture the optimal response, and this may depend on the type of input used. The abrupt changes in square wave inputs activate the transient overshoots (characteristic of signal termination in the autonomous system) in every input cycle (e.g., Fig. 5c, left), which are more prominent for the lower frequencies (compare Fig. 5b, left and right) and determine the output amplitude. These transients are not explicitly present in the responses to the gradual sinusoidal inputs. Studies with different input amplitudes and with the standard impedance (quotient of the Fourier transforms of the output and the input, in absolute value) are included in Supplementary Figures 4 and 5.

Fig. 5: Frequency-preference response in a CMC under periodic stimulation.

a Representative time courses for adapted sinusoidal inputs S_T (upper panels) with increasing frequencies and the corresponding outputs P (lower panels). For the three stimulation frequencies, the maximum reached by S_T is 35, whereas the minimum depends on the frequency, being 17.72, 24.81, and 30.01, in the left, middle, and right panels, respectively. The CMC parameter values are those listed for Fig. 4a and the periods of stimulation are, 33.33, 15.15, 6 s in left, middle, and right panels, respectively, resulting in frequencies of 0.03, 0.066, and 0.167 Hz. b Frequency response results, measuring amplitude (left) and gain (right) as a function of the input frequency. The amplitude of the output is defined as the difference between the maximum and the minimum values of P. The gain is the ratio between the amplitude of the output and that of the input. The three examples in a are indicated over the plot (same color code). c Representative time courses for the input S_T (upper panels) and the corresponding outputs P (lower panels), corresponding to an adapted train of square pulses stimulation. For the three stimulation frequencies, the maximum reached by S_T is 35, whereas the minimum depends on the frequency, being 20.3, 27.61, and 31.85 in the left, middle and right panels, respectively. Parameter values and input frequencies are as in a. d Frequency response results, measuring amplitude (left) and gain (right) as a function of the input frequency. The three examples in c are indicated over the plot (same color code).

We now turn to study in more detail the relationship between signal termination (in response to step inputs) and the preferred frequency responses of CMCs to oscillatory inputs. For each parameter set that exhibits signal termination (Fig. 3), we computed the output amplitude and gain profiles as described above (Fig. 5). We investigated the relationship between the two phenomena by comparing a number of representative attributes for the corresponding graphs. For the P vs t graph, we define t_dec and P_max (Fig. 6a, Constant Stimulation) as the time it takes P to decrease from its maximum to 63% of it and the maximum value of P, respectively. For the amplitude and gain profiles, we define (i) A₀ and G₀ as the amplitude and the gain obtained at the lowest frequency analyzed, (ii) f_A,max and A_max as the frequency and the amplitude at the preferred frequency (if there is one), (iii) f_G,max and G_max as the frequency and the gain at the preferred frequency (if there is one), and (iv) Q_A and Q_G as the ratios between the A_max and A₀, and G_max and G₀, respectively. In addition, we use the preferred periods T_A,max and T_G,max corresponding to the f_A,max and f_G,max (T_A,max = 1/f_A,max and T_G,max = 1/f_G,max). All these attributes are indicated in Fig. 6a and the definitions are included in the Methods section (Table 1).

Fig. 6: Parameter space exploration for periodic stimulation.

a Schematic representation of the two different scenarios analyzed, constant (left) and periodic (right) stimulation, and the attributes selected in each one. Constant stimulation is a step in S_T going from zero to the desired value at t = 0 and periodic stimulations are the adapted variations, sinusoidal, or train of square pulses, indicated in Figs. 4a–c and 4d–e, with frequencies in the range 10⁻⁴–10¹ Hz, resulting in periods in the range 0.1−10⁴ s. For constant stimulation, the attributes are t_dec and P_max. For periodic stimulations, the attributes are T_A,max, A_max, T_G,max, G_max (being T the period, T = 1/f). The frequency for which the maximum amplitude/gain is obtained and the corresponding maximum amplitude/gain is indicated over the plots as f_A,max, f_G,max, and A_max, G_max. b Results from the parameter sampling with 5,000 parameter sets. Sinusoidal stimulation was applied as S_T(t) = S_max(1+0.5sin(2πft))/1.5, then adapted as explained in Fig. 4 (and the text). Train of square pulses stimulation was applied as S_T = S_max for nT < t <(n+0.5)T, S_T = 0 for (n+0.5)T < t <(n + 1)T, n is a natural number, implying the input in the first half of the period. From the 5,000 parameter set values analyzed, 149 resulted in signal termination and all of them exhibited a maximum in the gain versus frequency plot, and also in the amplitude versus frequency plot but for sinusoidal stimulation only. Scatter plots relating an attribute from the output to periodic stimulation and an attribute from the output to step stimulation. The Pearson coefficient for each case is indicated over the plot. Left plots are T_G,max and T_A,max versus t_dec, right plots are G_max and A_max versus P_max. Scatter plots in dark yellow correspond to sinusoidal stimulation and scatter plots in cyan correspond to square stimulation.

Table 1 List and definitions of attributes analyzed in the article.

As for signal termination in Section 1, we define a number of criteria to establish the significance of preferred frequency response and we focus on the significant cases. First, both the amplitude and the gain are in the range 0–1 (the amplitude because of being normalized by S_T,max and the gain because P is part of the conservation condition of S_T stated in Fig. 1a, so P is always lower than S_T). We discard cases with a maximum amplitude/gain lower than 0.1. Second, we select those outputs for which A₀ < 0.9 A_max and G₀ < 0.9 G_max. This ensures a significant gain with a large enough maximum. To characterize the strength of the preferred frequency response, we define the peak-to-initial value Q_A = A_max/A₀ in the amplitude profiles, and Q_G = G_max/G₀ in the gain profiles. This implies Q_A, Q_G > 1.1. While the specific numbers chosen (0.1, 0.9) are somehow arbitrary, similar results to the ones presented in this paper are obtained for other choices. We refer to the parameter sets leading to these two conditions under periodic stimulation, as having frequency preference.

The results of the parameter space exploration indicate that for all the parameter sets for which the CMCs show signal termination, they also show a preferred frequency response in both the amplitude and gain to sinusoidal stimulation. For square wave stimulation, this conclusion holds only for the gain since the amplitude does not show any preferred frequency response. Furthermore, we observed no cases where CMCs exhibit preferred frequency responses without signal termination. In Fig. 6b, we present the correlation graphs between the attributes for signal termination (abscissas) and the preferred frequency responses (ordinates). We include those exhibiting higher correlations, as measured by the Pearson coefficient. The correlation between the decay time and the preferred periods (Fig. 6b, left) indicates that the effective slower time scales of the autonomous CMCs play a significant role in determining their preferred responses to sinusoidal inputs. Unlike linear systems, these time scales are not explicit as model time constants, but involve a combination of model parameters. These results show that input is more effective in producing a significant output if it allows the system to terminate the signal before the new stimulation cycles begin. The correlations between A_max and G_max confirm the relationship between the time responses of CMC cycles to constant pulses and the frequency responses to oscillatory inputs. Although these results may not be completely unexpected, they are also not obvious for CMCs, particularly taking into account that CMCs are not explicitly designed as negative feedback loops.

We now turn to the analysis of the biochemical mechanisms of generation of preferred frequency responses in the gain and we compare them with the mechanisms of signal termination to further strengthen the relationship between the two phenomena. These mechanisms consist of the transition from a low-pass filter to a band-pass filter as the parameter under consideration increases. We have identified a number of paths in parameter space that lead to the generation of preferred frequency responses in the gain. They are summarized in Fig. 7. They involve the changes in the quotient of two model parameters (all other parameters fixed) and therefore each of these paths is degenerate in the sense that multiple combinations of model parameters satisfy the same constant constraint for the quotient (Fig. 7b, c).

Fig. 7: Analysis of the mechanisms leading to frequency preference.

a P versus time for step (constant) stimulation. b Gain versus frequency under sinusoidal stimulation. c Gain versus frequency under the train of square pulses stimulation. Each row corresponds to different parameter variations controlling the transition in or out of the signal termination region. First row, E_1T/S_T; second row, E_2T/S_T; third row, k₁ = αk_2ref; fourth row, k₂ = βk_1ref, with k_2ref and k_1ref the corresponding parameter values for a reference set (as in Fig. 5). These parameter variations are color-coded as indicated in the color bar at the bottom of the figure.

In principle, from a low-pass filter, preferred frequency responses in the gain profile can be generated as the result of the increase of a parameter value (or combination of parameter values) because (M1) G increases faster for intermediate frequencies than for lower frequencies, or (M2) G decreases faster for lower frequencies than for intermediate frequencies.

Increasing values of E_1T/S_T generate preferred gain responses by mechanism M1 (Fig. 7b, c, row 1) and it is more pronounced and robust for sinusoidal than for square wave inputs (compare panels b and c). For the latter, there is a transition to a low-pass filter as E_1T/S_1T continues to increase.

In contrast, increasing values of E_2T/S_T generate preferred gain responses by mechanism M2 (Figs. 7b, c, row 2). As for the previous case, the resonances are more pronounced and robust for sinusoidal than for square wave inputs (compare panels b and c), and for the latter, there is a transition to a low-pass filter as E_2T/S_2T continues to increase. Therefore, decreasing values of E_2T/S_T also generate preferred frequency responses in gain for square wave inputs by mechanism M1.

Increasing values of α generate preferred gain responses by mechanism M1 (Fig. 7b and c, row 3), which are also more pronounced and robust for sinusoidal than for square wave inputs (compare panels b and c) and, for square wave inputs the preferred frequency response is terminated as α continues to increase. Finally, decreasing values of β generate preferred gain responses by mechanism M2 (Fig. 7b and c, row 4), which are also more pronounced and robust for sinusoidal than for square wave inputs (compare panels b and c).

Cascades of CMCs

In this section, we consider a cascade composed of two CMCs, as indicated in Fig. 8a where the first cycle (CMC₁) is subject to constant or periodic stimulation and we examine how these signals are processed across the network. More specifically, we first study whether signal termination (under constant stimulation) is affected by the presence of a downstream cycle, if it can be propagated downstream in the cascade, and if new behavior emerges from the coupling of cycles. We then characterize the frequency response of the cascade under periodic stimulation. For simplicity, we restrict our study to the gain in response to sinusoidal stimulation and we leave out the details corresponding to the amplitude. We emphasize that a cascade is not simply a feedforward network where CMC₂ responds to the output from CMC₁, but it is an interconnected network where CMC₂ affects CMC₁ and therefore knowledge from the CMC₁ output is not enough to predict the CMC₂ output.

Fig. 8: Different regimes in cascades of CMCs.

a Scheme representing a cascade with two CMC cycles. S₁ and S₂ are the unmodified substrates, and P₁ and P₂ are the products (modified substrates), in each level. E₁ and E₁’ are the modifying enzymes or kinase and phosphatase in the upper level, C₁ and C₁’ are the complexes substrate–kinase and substrate-phosphatase in the upper level, P₁ and E₂’ are the modifying enzymes in the lower level, C₂ and C₂’ are the complexes in the lower level. The three kinetic parameters associated with each enzymatic reaction are a_i, d_i, k_i, the association, dissociation, and catalytic rates, respectively (more details in the Methods section). The conservation conditions are listed in the scheme. b Table summarizing the five representative numerical studies carried out. In all of them, step stimulation is applied. In studies 1 and 2, parameter values of the first cycle are selected so that this cycle in isolation exhibits signal termination or monotonic response, respectively, whereas the parameter values of the second cycle are randomly selected from predefined ranges. In studies 3 and 4, parameter values of the second cycle are selected so that this cycle in isolation exhibits signal termination or monotonic response, respectively, whereas the parameter values of the first cycle are randomly selected from predefined ranges. In Study 5, all the parameter values, first and second cycle, are randomly selected from the mentioned ranges. In all, 5,000 results of numerical simulations were analyzed in studies 1–4, and 10,000 in Study 5. Parameter values for each study are listed in Supplementary Table 1. The ranges for random exploration are concentrations: 10–100, association and dissociation rates: 0.1–10, and catalytic rates: 1–50. The results were classified into four groups (right): (i) those exhibiting signal termination in the first cycle but not in the second one (red bar); (ii) the inverse of (i) (light blue bar); (iii) those exhibiting signal termination in both cycles (orange bar); (iv) those with monotonic behavior in both cycles (light green bar). These groups are schematized next to the table (with the same color pattern as in the last row of the table). The criteria used to classify the outputs in groups are described in Methods. The lowest row in the table contains the results of each study, with a bar plot representing the percentage of cases in groups i–iv.

We summarize our studies and results for signal termination in Fig. 8b. In the first two studies (study 1 and study 2), we selected specific parameter values for the first cycle (CMC₁) and the parameter values for the second cycle (CMC₂) were randomly distributed. In the following two studies (study 3 and study 4) we selected specific parameter values for the second cycle (CMC₂) and the parameter values for the first cycle (CMC₁) were randomly distributed. In study 1 or 3, CMC₁ or CMC₂, respectively, exhibits signal termination (under constant stimulation) and preferred frequency response (periodic stimulation). In study 2 or 4, CMC₁ or CMC₂, respectively, exhibits no signal termination and therefore no preferred frequency response. In study 5, the parameter values of CMC₁ and CMC₂ were randomly selected. There are four possible types of responses for the cascade (Fig. 8b, right). The last row in the table indicates the percentage of each of them. We briefly report our results below and present a comprehensive table comparing the results of the five studies in Supplementary Figure 11).

In study 1, we found that ~83% of the explored parameter sets (initial exploration: 5,000 sets, significant responses: 4,212, see Methods for a detailed description of the criteria) resulted in the persistence of signal termination for CMC₁, whereas CMC₂ responds monotonically. A smaller number of cases resulted in signal termination in both cycles (~15%). These numbers indicate that signal termination is eliminated ~2% of the time. Why can the second cycle reduce the signal termination property in the first cycle?. The presence of the second cycle introduces a competitive inhibition effect for the reverse reaction in the first cycle, i.e, the one converting the product back to substrate. This competitive inhibition results in an effective reaction velocity with an increased Michaelis–Menten constant Km³¹. Since the phosphatase affinity for the product is the inverse of this Km, the presence of the second cycle decreases the phosphatase affinity, whereas the kinase affinity stays unmodified. As indicated in Fig. 3, a high ratio Aff₂/Aff₁ promotes signal termination, so the competitive inhibition effect that effectively reduces Aff₂ results in a decrease in the number of cases with signal termination.

In study 1, most but not all of the cases with signal termination in CMC₁ resulted in frequency preference when periodically stimulated (65%). This is different from our expectation from the responses of individual CMCs to periodic stimulation (where all cases showing signal termination also showed frequency-preference response). Together, although both signal termination and frequency preference are not always propagated by the cascade, they persist in a non-negligible number of cases.

In study 2, we found that the absence of signal termination in the isolated CMC₁ does not necessarily cause the cascade to show a monotonic response. In fact, we found signal termination in CMC₁ only in 26% of the cases. In this case, sequestration of the first cycle substrate in its modified form, i.e., the product, by the substrate of the second cycle reduces the available first cycle substrate for the first cycle kinase and phosphatase, leading to an enzyme-in-excess regime. This situation is not general and depends on the parameter values, resulting in 26% of signal-termination new cases. An alternative explanation is that by adding the second cycle there is an extra term of sequestration of the product of the first cycle, contributing in this way to terminate the signal (see Supplementary Figure 10). Although signal termination in CMC₂ is not necessarily unexpected taking into account that a monotonic response can rapidly reach steady-state, therefore mimicking an effective constant stimulation situation to CMC₂, particularly when the isolated CMC₂ shows signal termination, signal termination in CMC₁ is a network effect.

In study 3, we found that adding an upstream cycle does not remove signal termination in CMC₂. Study 4 indicates that is very unlikely that an upstream cycle can lead CMC₂ into signal termination when isolated produces a monotonic response. In study 5, we also found signal termination and frequency preference in both cycles (signal termination: 14% in CMC₁, 4% in CMC₂, 0.7% in both; frequency preference: 5% in CMC₁, 0.6% in CMC₂, 0.44% in both).

Of particular importance are the network processing capabilities of cascades in controlling the properties of signal termination and frequency-preference responses (to constant and oscillatory inputs, respectively). To illustrate these issues, here we focus on those outputs from study 1 exhibiting signal termination and frequency preference in both cycles of the cascade. In Fig. 9, we present examples of both cascade amplification and attenuation in both signal termination (Fig. 9a) and frequency preference (Fig. 9b). We used the following metrics for each CMC: Q = P_max/P_ss and Q_G = G_max/G_0.

Fig. 9: Amplification and attenuation in cascades of CMCs.

a Representative example of time courses exhibiting amplification (left) or attenuation (right) when stimulated with a step. b Representative examples of gain versus frequency plots exhibiting amplification (left) or attenuation (right) when periodically stimulated.

Signal termination and frequency-preference response is not captured by usual CMC’s approximations

Modeling studies of enzymatic reactions often rely on a number of simplifying assumptions. Perhaps, the best known is the Michaelis–Menten approximation, which reduces enzymatic reaction models to single differential equations (one-dimensional systems). A key assumption is the identification of a fast time scale governing the evolution of the initial, transient increase in the complex concentration, which then relaxes to zero on a much slower time scale giving rise to the product increase^35,36. This dimensionality reduction approximation affects minimally the dynamics of the product increase provided a good estimation of the complex concentration after the short transient. For this and other systems where the dynamics are quasi-one-dimensional, the elimination of the transient dynamics is practically inconsequential. In contrast, their response to external inputs fails to capture the response of original systems. Two prototypical examples are signal termination (in response to constant inputs) and the preferred frequency responses to oscillatory inputs. In linear systems, for example, these two phenomena are present in two-dimensional models for the appropriate parameter regimes (e.g., slow negative feedback effect), but they are absent for their quasi-one-dimensional approximation where the autonomous transients are neglected. For signal termination, the explanation is rather simple. For signal termination to occur, the temporal profile of the measured variable achieves a maximum prior to reaching the steady-state (and different from it). For this to happen, the action of a second variable opposing the increase is necessary. Otherwise, the principle of uniqueness for well enough behaved (one-dimensional) systems would be violated, and therefore the steady-state and maximum coincide. For the preferred frequency response, the necessity of higher-than-one dimensionality is predicted by standard calculations. The dynamic explanation is more involved and it derives from the dynamical system analysis carried out^24,37. These two phenomena are not observed in single enzymatic reactions, as they lack the main ingredients (slow negative feedback), but, as we showed, they are present in networks of enzymatic reactions.

Our study in this paper involves the mechanistic formulation of CMCs (and cascades) without any simplifying assumption. In this section, we examine the effectiveness of the simplifying assumptions typically used for modeling CMCs in reproducing signal termination and frequency preference. We do not intend to carry out a systematic analysis, but rather to provide some insight into the issues and the motivation for our modeling decisions.

The QSSA is also frequently used to derive reduced models for enzyme-catalyzed reaction networks (as CMCs). Under the additional substrate-in-excess assumption (the substrate is in excess over the enzymes), Goldbeter and Koshland derived a single differential equation for the modified substrate P¹ (details in the Supplemental Information). The applicability of this last model is restricted to conditions when the substrate concentration is much higher than that of the converter enzymes. In order to correct for the failure of the approximation to capture the transient dynamics, another approximation called the tQSSA was designed³⁸. It is based on certain linear combinations of the original variables and it has proven to yield much better approximations, especially when the enzyme concentrations become comparable to or larger than that of the substrate⁶. However, the description of a CMC with a tQSSA also leads to a single differential equation for the modified substrate P (details in the Supplemental Section 7), and therefore it is unable to capture both signal termination and the frequency-preference response to oscillatory inputs (Fig. 10).

Fig. 10: Approximation comparisons.

a Time courses for variable P obtained by using the mechanistic model (black solid), the tQSSA (orange dashed), and a combined mechanistic model with a Goldbeter–Koshland approximation for only one of the modified enzymes (purple and green). b Gain versus frequency plots for the same conditions as in a. Parameter values as in Fig. 5.

To examine whether a mixed type of approximation that reduces the model dimensionality, but results in a two-dimensional (rather than one-dimensional) model, is able to capture both signal termination and the frequency preferred responses, we combined the mechanistic model with a Goldbeter–Koshland type of approximation for only one of the enzymes leading to two reduced models (C₁ GK and C₂ GK), one of each enzyme (E₁ and E₂, respectively). Our results (Fig. 10, purple and green lines) demonstrate that while these reduced models show signal termination and frequency-preference responses, the resulting patterns are not good approximations for the full, mechanistic model.

Discussion

CMCs and cascades of CMCs have been studied both theoretically and experimentally as the primary intracellular signaling mechanisms in living systems^{1,3,7,17,18,29,39,40}. An important question is how these systems respond to external stimulation. Two prototypical stimulation signals are step-constant (abrupt transitions between zero stimulation and a constant positive stimulation) and periodic. Step-constant stimulation serves the dual purpose of characterizing the variety of steady states available to CMCs and uncovering the transient dynamics leading to these states. In the simplest case, these dynamics are monotonic. However, recent work has shown that CMCs have the ability to exhibit signal termination¹⁴ (or its counterpart, adaptation) where the initial response of the system is higher than the steady-state level, or, from another perspective, the system has the ability to produce a stronger response than that predicted by the steady-state, but only for a restricted amount of time. This effect is analogous to what is called adaptation in sensory and biochemical systems, and resembles an overshot in the temporal response.

Signal termination is a relatively recent finding in CMCs, which has been observed in the enzyme-in-excess regime, but not in the substrate-in-excess regime, which has been the most studied so far. Periodic stimulation of CMCs serves the purpose of understanding how these systems react to time-varying inputs. On the one hand, periodic stimulation is arguably the simplest type of non-decaying (with time) input. On the other hand, the response to sinusoidal inputs can be used to reconstruct the response signal to a rather general class of time-varying inputs. Square wave inputs constitute a link between both the step-constant and sinusoidal inputs where the input is periodic and changes abruptly between two constant values, thus periodically reproducing the transient effects described above in a frequency-dependent manner. The response of CMCs to periodic stimulation has been primarily studied in the substrate-in-excess regime where the product responses have been shown to be low-pass filters. However, outputs that terminate in response to step-constant inputs have been linked to band-pass filter responses to sinusoidal inputs^24,25,26,27, thus raising the possibility that CMCs in the enzyme-in-excess regime exhibit this type of responses to periodic inputs in general, and both signal termination and preferred frequency responses are communicated across CMCs in the cascades of which they are part. In this paper, we set out to address these issues via a combined modeling and computational approach.

It is common practice to consider that the number of enzymes, either kinase, or the ratio between kinase and phosphatase, drive the activity in a CMC^3,10,28,41. This is typically done in the substrate-in-excess regime. However, in this article, we considered the opposite, enzyme-in-excess regime. In this case, stimulating the enzyme levels produces almost no effect as shown in Supplementary Figure 1). We then stimulated the system by changing the amount of substrate, particularly, the total level of the substrate. It is also important to notice that in several cases substrates are enzymes as well. This is the case of kinase cascades in which the activated substrate in one level works as the enzyme for the downstream level. The MAPK cascade is a well-known example in this regard¹⁰. In these cascades, stimulating the enzymes or the substrates could play interchangeable roles.

The total amount of proteins in a CMC usually changes in a slower time scale than that of signaling, but that is not always the case and several studies consider total variations in the same time scale as that of signaling^29,39. The effect of regulating the total level of a signaling protein and the impact this could have in its dynamics was comprehensively studied both theoretically and computationally^22,42. Even when experimental results are available in the analogous ligand–receptor reaction scenario, the common practice is to consider that the ligand is the stimulus. However, in certain articles, it was convenient to stimulate the receptor level⁴³.

In the same way, the total amount of substrate was treated as a signal in this paper, other factors can be treated as a signal too. In Supplementary Figure 2 the signal is k₁, the kinase catalytic rate. A step-like variation in k₁ leads to a signal termination response, while a periodic forcing in k₁ leads to a frequency-preference response. Similar studies are found in the literature considering different factors as a signal, as the substrate unbinding rate⁴⁴, the kinase inactivation rate⁴⁵, and the degradation rate⁴⁶, to mention some examples.

The stimulation selected in this study, the total substrate concentration, places this article among those evaluating absolute concentration robustness (ACR)^47,48. ACR implies that the steady-state output of a system is completely independent of a protein concentration. For CMCs with high enzyme concentration, i.e., CMCs in the enzyme-in-excess regime, it was reported that concentration robustness only occurs under certain conditions with respect to total enzyme concentrations⁶. Even though we have not evaluated theoretically ACR in the enzyme-in-excess regime with respect to total substrate concentration, the results we obtained indicate that it does not occur. The reason for this is that systems with ACR have robust perfect adaptation for the target species. A system possesses perfect adaptation if the output of the system returns to the pre-stimulus level after the value of the input parameter is changed to a different constant level^{49,50,51,52,53}. If the system achieves perfect adaptation independently of the system parameters, it is said to have robust perfect adaptation. As shown in this article, after the transient response the system does not generally return to pre-stimulation levels. In summary, CMCs in the enzyme-in-excess regime exhibit neither perfect adaptation nor ACR with respect to total substrate concentration in a purely kinetic framework. However, it could be possible that when considering spatial organization, even though the total enzyme may be in excess, the dynamics might be different locally⁵⁴.

Motivated by the predictions in¹⁴, we characterized the signal termination responses in terms of a number of attributes (peak P_max and steady-state P_ss) and excluded from the analysis those cases for which signal termination is mild (small P_max, low P_max/P_ss). The key elements for the signal to terminate are the relative amounts of enzymes and substrate, the relative velocities of the two modifying enzymes, and their relative affinities. The underlying mechanism is relatively simple, when the substrate is added in a step-like profile, a fast kinase modifies it producing the product, while a slower phosphatase undergoes the opposite process terminating the signal. Our results (Figs. 1–3) quantitatively confirmed the conditions stated in (Bluthgen et al., 2006), namely that if the enzymes are in large excess over the substrate, if the kinase is fast and has low affinity, and if the phosphatase is slow and has high affinity, then signal termination occurs. We referred to them as the Bluthgen conditions.

To more deeply understand and further characterize signal termination we combined analytical calculations, local analysis, and a parameter space exploration. Some of the studies were performed over specific, representative parameter sets (Figs. 1, 5, and 7), whereas others are based on this parameter space exploration (Figs. 3, 6, and 8). The representative parameter sets capture scenarios that have more general validity. For this first group, we included the results of alternative parameter sets in the Supplementary Information, showing that the conclusions extracted from those studies are independent of the choice of parameter set (Supplementary Section 5).

The goals of the parameter space exploration were: (i) to find out the connection between parameter sets satisfying Bluthgen conditions and outputs exhibiting signal termination (according to the criteria defined in Section 1); (ii) to investigate whether signal termination persists if any of the four Bluthgen conditions (E_1T/S_T > 1, E_2T/S_T > 2, V₁/V₂ > 1, Aff₂/Aff₁ > 1) are relaxed, which one can be relaxed, to what extent, and which exert tighter control over signal termination; and iii) to understand which parameter set characteristics lead to strong signal termination. Strong signal termination is a more restrictive condition in which the signal terminates to almost pre-stimulation levels, close to what is known as perfect adaptation.

Although signal termination arises in an enzyme-in-excess regime, our results show that this is not a necessary condition, but signal termination can occur in a regime where the enzymes and the substrate are in similar (comparable) amounts. Particularly, the phosphatase level has a tighter control on signal termination, being this enzyme is the one that converts back the product into the unmodified substrate thus terminating the signal. The velocities and affinities requirements (V₁ > V₂ and Aff₂ > Aff₁) can be relaxed while still observing signal termination, but one at a time. In fact, relaxing both together leads to a monotonic temporal profile. If condition E_2T/S_T > 1 is not satisfied, conditions Aff₂/Aff₁ and V₁/V₂ are both required, instead, if it is condition E_1T/S_T > 1 that is not satisfied, only one extra condition (Aff₂/Aff₁ or V₁/V₂) can lead to signal termination (this conclusion is extracted from Fig. 3c).

Before we could determine whether the parameter sets leading to signal termination under constant stimulation ensured a preferred frequency response under periodic stimulation, we needed to overcome an important issue, which is common to any system in which the quantity that needs to be controlled satisfies a conservation law and is under sequestration conditions: S_T cannot follow the desired periodic profile, an adapted but still, periodic one was used instead (Fig. 4). This must be done so as not to include negative values for amounts of the participating species. Specifically, the chemical species at play cannot follow any desired periodic profile, the time scale involved in the release from sequestration together with the forcing time scale determine an adapted stimulation profile. This led to a modification of the stimulation protocol that sacrificed uniformity across frequencies in the input for the sake of understanding the preferred frequency properties of CMCs under realistic conditions and appropriate normalization that made the comparisons across frequencies possible.

The adaptation of the input signal is a necessity dictated by the system and not our choice. Our choices lie simply in (1) the type of adaptation among more than one possible option, and (2) the choice of a system that presents this problem. One solution to the latter would have been to avoid the study of this system and choose a “more elegant one”. However, we believe we are shedding light on a problem that presents itself in realistic situations. To our knowledge, frequency-preference studies of complex enzymatic reactions are not abundant in the literature, most likely because of the type of issues we are discussing here. In the cases where results were successfully obtained without the need for adaptation in the input signal, it is not clear whether all the variables remain positive or only the variables presented in the graphs and whether the model reductions used conserve all the frequency-dependent information present in the full model from where they originate. To our knowledge, these issues have not been discussed in the literature on CMCs and related biochemical systems, and call for further research and examination.

The main idea we test in this paper is whether there is a relationship between signal termination and the amplitude profile properties of CMCs to periodic inputs. This knowledge allows inferring the presence of time scales that support these phenomena. While engineering studies of linear systems focus on prediction, other scientific studies (e.g., cell signaling and neuroscience) focus on inferring the presence of time scales that are common to different, apparently disconnected phenomena. In this paper, we follow this tradition. We acknowledged and discussed in detail the difficulties in applying a uniform protocol across frequencies. However, by normalizing by the input amplitude we overcome this issue, if not totally, partially, but this does not affect our ability to learn from the system what we intended to learn. Furthermore, Supplementary Figure 4 shows that the gain is not affected by the amplitude of the input, a result that was expected by the definition of gain. In addition, the responses measured in terms of the standard impedance (quotient of the Fourier transforms of the output and the input, in absolute value) have predictive value, as indicated in Supplementary Figure 5.

Regarding which biological situation could lead to such a periodic control of the total amount of substrate, it is well-known that signaling proteins are promiscuous, i.e., the same protein is shared by different pathways leading in several cases to a specificity problem^55,56. The use of the free amount of protein by others pathways, particularly oscillating pathways could lead to the situation described in this article. Although we agree, these pathways will most likely not lead to nice sinusoidal or square wave inputs, these are, in our opinion, the most sensible ones for an organized study of the phenomena.

We evaluated the existence of preferred (optimal) frequency responses to periodic stimulation (band-pass filters) in terms of two quantities that provide complementary information: amplitude and gain. We used the two types of input protocols mentioned above, which also provide complementary information: sinusoidal and square waves. We found that the CMC exhibited preferred frequency responses in the gain for parameter values for which the underlying autonomous system exhibits signal termination, whereas the output amplitude did not necessarily capture the optimal response, and this depended on the type of input protocol. Analyzing correlation graphs between the attributes for signal termination and corresponding preferred frequency responses, we found that the effective slower time scales of the autonomous CMCs played a significant role in determining their preferred responses to sinusoidal inputs. These results indicated that input is more effective in producing a significant output if it allows the system to terminate the signal before the new stimulation cycle begins. These results are consistent with previous findings in neuronal systems²⁴.

From the biochemical point of view, the preferred frequency response in the gain profile can be generated by two different mechanisms as the relevant biochemical parameters change (Fig. 7): (i) G (gain) increases faster for intermediate frequencies than for lower frequencies, and (ii) G decreases faster for lower frequencies than for intermediate frequencies. These results together with the results discussed above strengthen the relationship between the emergence of preferred frequency response to periodic inputs and signal termination beyond the parameter exploration exercise.

Testing the hypothesis that signal termination and preferred frequency response to oscillatory inputs in isolated CMCs play a role for the responses of the networks in which these cycles are embedded, to both step-constant and periodic stimulation, requires the identification of representative case studies to avoid the rapid explosion of parameter combinations. In this paper, we focused on the minimal network model consisting of a cascade of two CMCs and identified five cases studies that allowed us to evaluate whether the signal is affected by a downstream cycle and if it can be propagated downstream in the cascade (Fig. 8, Study 1), if this behavior can emerge from the coupling of cycles as a network effect (Fig. 8, Studies 2–5), and if new behavior can emerge from the coupling of CMCs. Within the explored ranges we found no new dynamic behavior, contrary to our expectation of uncovering damped oscillations⁷. We concluded that signal termination can be propagated to a downstream cycle and that this second cycle is not able to remove the signal termination in the upper one.

Interestingly, we found that a monotonic behavior in the first cycle can be transformed into signal termination by the addition of a second cycle with particular parameter values (Fig. 8, Study 2, output i). The opposite study, i.e., sampling on the upper cycle, is not able to remove signal termination in the downstream one (Fig. 8, Study 3), reinforcing that the control of the behavior is exerted by the phosphatase of the cycle with signal termination, cycle 2 in the considered study. These studies on cascades were complemented with the search of preferred frequency responses and the analysis of amplification and attenuation in the cascade (Fig. 9).

The nature of the perturbation of the first cycle by the second that alters the first cycle’s behavior is called retroactivity, i.e, the modulation of the first cycle’s behavior performed by its load, the second cycle^7,17,18,57. This change of behavior can have different effects depending on the parameter conditions.

The discussion on simplifying assumptions for modeling multi-step pathways in systems biology is very active in the literature⁵⁸. Based on it, one would wonder if it is possible to attempt some theoretical analysis of the response of the system. As we see it, this type of analysis could go in three directions: (i) reducing the dimensionality of the system to two dimensions, (ii) linearizing the system and preserving the three dimensions, or (iii) developing tools to analyze the original 3D system. The reduction of dimensionality necessarily needs to be based on the identification of a separation of time scales and keeping only the portion with slow dynamics. But, we in fact show in the paper that the relevant dynamics are lost if we proceed in this way (Section 4, signal termination and frequency-preference response are not captured by usual CMC’s approximations). Linearizing the system around the equilibrium would produce the expected results about signal termination and resonance, but its validity will be restricted to a very small region in the phase-space around this equilibrium, and therefore it would not be a reasonable analytical approach to explain our results. The development of analytical tools to analyze the original 3D system is the most intellectually tempting option and there are various possible approaches. One of them is using piecewise linear systems and analyzing the evolution of the trajectory as it crosses regimes with different combinations of linear pieces. Another one is to track the responses of (“moving”) 3D phase-space diagrams in response to the inputs with the aid of the graphs of null surfaces. Both of these approaches are complex and are beyond the scope of the present paper.

An important feature of our study is that our (mechanistic) modeling approach did not use the available and ubiquitous approximations, which often render one-dimensional models that are expected to be unable to produce signal termination and preferred frequency responses to oscillatory inputs. To confirm this and to further understand the effects of the removal of the fast time scales from the models by the application of the dimensionality reduction processes, we repeated some of the protocols using the reduced approximated models. As expected, the tQSSA produced no signal termination and low-pass filtering. The other approximations considered did produce signal termination and preferred frequency response to periodic inputs, but not good approximations to similar results using the mechanistic model (Fig. 10). Analogous consequences are found in other areas, in which the search of simple and compact models leads to dimensionality reduction that precludes non-monotonic temporal profiles such as the ones described in this article^59,60.

Different network topologies can lead to an overshoot transient response. The incoherent feedforward loop is one of them¹⁹. However, the ability of a CMC to produce similar transient responses was only studied by Bluthgen and colleagues. Furthermore, the groundbreaking article by Ma and colleagues defining network topologies that can achieve biochemical adaptation¹⁹, considers three-node topologies where each node is a CMC. This indicates that the literature considers much more complicated systems than an isolated CMC in order to explain overshoot/signal termination/adaptation responses. The unstated implication is that these topologies are the minimal ones that produce the phenomenon. We reasoned that if the simple two-dimensional linear systems are able to show signal termination, then CMCs should also be able to display this phenomenon, at least in a regime where the dynamics unfold in more than one effective dimension and the necessary time scales are preserved. In fact, our studies indicate that the reason behind the current understanding in the literature is twofold: (1) CMCs were mostly studied in the substrate-in-excess regime, where signal termination does not occur, (2) signal termination is not captured by the usual approximations, as they rely on the removal of the fast time scales. These observations partly explain why signal termination and frequency-preference response were overlooked in this ubiquitous signaling motif.

Regarding the modeling approach in this article, there are instances where modification cycle/cascades may not be described in ODE terms and these include various contexts where the spatial dimension is important⁶¹. This issue requires more research.

Summarizing, substrate sequestration by its modifying enzymes, and in particular by the phosphatase, might be a means to achieve signal termination and desensitization downstream of receptors without involving an explicit negative feedback loop. This behavior can be at play under physiological conditions, as the abundances of substrates and enzymes are similar in vivo. Therefore, it is relevant to study how CMCs, cascades of CMCs, and signaling pathways combining them, process different temporal inputs under sequestration conditions. Our predictions about the band-pass filter behavior of these elements (when periodically stimulated) shed interesting insights into fundamental biological processes and the computations that might be carried in biochemical networks. More research is needed to test these predictions experimentally and to examine the functional consequences of band-pass filter behavior in CMCs for larger networks in which CMCs are embedded. Although preferred frequency responses of CMCs could be an epiphenomenon, research in other areas shows they may be functional for the generation of network oscillations^62,63.

Methods

All the algorithms used in this paper were written in Python 3.8, Spyder 4. The libraries used were numpy, numba, pyDOE, pylab.

Mechanistic description for a CMC

Following the scheme, reactions, variables, and parameters indicated in Fig. 1a and using the Law of Mass Action, the ode system describing a CMC is the following:

$$\frac{{\rm{d}}{C_1}}{\rm{dt}} = a_1E_1S - \left( {d_1 + k_1} \right)C_1$$

$$\frac{{\rm{d}}{C_2}}{\rm{dt}} = a_2E_2P - \left( {d_2 + k_2} \right)C_2$$

$$\frac{{\rm{d}}{P}}{\rm{dt}}= k_1C_1 - a_2PE_2 + d_2C_2,$$

combined with the following conservations laws:

$$E_{1T} = E_1 + C_1;\;E_{2T} = E_2 + C_2;\;S_T = P + S + C_1 + C_2$$

Mechanistic description for a cascade with two CMCs

Following the scheme, reactions, variables, and parameters indicated in Fig. 8a and using the Law of Mass Action, the ode system describing a two-level cascade is the following:

$$\frac{{\rm{d}}{C_1}}{\rm{dt}} = a_1E_1S_1 - \left( {d_1 + k_1} \right)C_1$$

$$\frac{{\rm{d}}{C_1^\prime}}{\rm{dt}}= a_1^\prime E_1^\prime P_1 - \left( {d_1^\prime + k_1^\prime } \right)C_1^\prime$$

$$\frac{{\rm{d}}{C_2}}{\rm{dt}}= a_2P_1S_2 - \left( {k_2 + d_2} \right)C_2$$

$$\frac{{\rm{d}}{C_2^\prime}}{\rm{dt}}= a_2^\prime E_2^\prime P_2 - \left( {k_2^\prime + d_2^\prime } \right)C_2^\prime$$

$$\frac{{\rm{d}}{P_1}}{\rm{dt}}= k_1C_1 - a_1^\prime P_1E_1^\prime + d_1^\prime C_1^\prime - a_2P_1S_2 + \left( {d_2 + k_2} \right)C_2$$

$$\frac{{\rm{d}}{P_2}}{\rm{dt}}= k_2C_2 - a_2^\prime E_2^\prime P_2 + d_2^\prime C_2^\prime$$

combined with the following conservations laws:

$$E_{1T} = E_1 + C_1;\;E_{1T}^\prime = E_1^\prime + C_1^\prime ;\;E_{2T}^\prime = E_2^\prime + C_2^\prime$$

$$S_{1T} = S_1 + P_1 + C_1 + C_1^\prime + C_2;\;S_{2T} = S_2 + P_2 + C_2 + C_2^\prime$$

Parameters definitions and units

a, association rate

d, dissociation rate

k, catalytic rate

E_T, total amount of kinase

E_T’, total amount of phosphatase

S_T, total substrate

Rates a, d, k could appear without or with an apostrophe, indicating that they are associated with the reaction catalyzed by the kinase or the phosphatase, respectively. Those rates could appear with a subindex 1 or 2, indicating the first and second cycle in the case of a cascade. For the numerical simulations, the initial concentrations are such that all the substrate is in the unmodified form and all the enzymes are in the free form (i.e., are not forming complexes). The values selected for the parameters are indicated in each figure and summarized in Supplementary Table 1.

Some of the parameters have dimensions of concentration (E_T, E_T’, S_T), some have dimensions of 1/time (d, k), and others have dimensions of 1/(time · concentration) (a). The unit of time is selected as minute (so that the duration of the responses in signal termination is consistent with significant biological durations), and the unit of concentration is arbitrary. This means that when a value for a parameter is listed, the corresponding unit has to be added, for example, if the value of parameter a is 10, it has to be read as 10 × 1/(min × concentration). Once the reference concentration is selected, all the parameters follow that reference unit. The interpretation of the results depends on the choice of the reference unit concentration (the ‘0’ in log scale). For example, if the reference dimensional concentration is chosen as 0.1 μM, this leads to interpreting the scanned intervals (from 10¹ to 10²) as being in the range (1 μM, 10 μM). However, this is just an example and the choice of the reference unit concentration remains a degree of freedom in our numerical methodology.

Parameter space exploration for signal termination

We took as an initial guide for parameter values and ranges, the intervals obtained from the literature and used in previous articles^40,64 and the table published in the Supplementary information of Bluthgen’s article¹⁴. The first intervals indicated a range of 0.3–224 for the enzymes and 3–1000 for the substrate, evidencing that they were based on a substrate-in-excess regime. We modified them for this article so that both enzymes and the substrate take values in the same range, 10–100. k₁, k₂ are in the range 1–50, which is similar to that in the mentioned intervals (6.3–600). Parameters a₁, d₁, a₂, d₂ were selected in the range 0.1–10 according to Bluthgen’s table. Every parameter set in the paper was taken from these ranges.

We logarithmically sampled parameters a₁, d₁, a₂, d₂ between 10⁻¹ and 10¹, k₁, k₂ between 1 and 50, and S_T, E_1T, E_2T between 10¹ and 10² using Latin Hypercube Sampling⁶⁵ for 20,000 different sets. To do that, we used pyDOE Python library and sampled linearly in the exponents. We scanned the 20,000 sets, dividing the process into four independent routines of 5,000 sets of different samples.

Parameter space exploration for frequency preference

The frequency response was studied in 5,000 parameters sets. Each set was stimulated with 10 different frequencies f, equally spaced logarithmically in the range 10^–2/60–10–1/60, with step 10^0.1/60. In each case, gain and amplitude were measured. If the profile obtained with these ten data points indicates an increasing function, it is concluded that it will exhibit a frequency preference, so f is increased until finding a decreasing trend. In this process, the preferred frequency is obtained. If instead, the profile (gain or amplitude versus frequency) shows a decreasing trend for all the explored frequencies, it is assumed that this is a low-pass filter. G₀ always corresponds to G₀=G(f = 10⁻²/60). Once a frequency-preference case is identified, a finer grid is used to better capture f_max (f is varied between 0.5f_max and 2f_max with a step of 0.05f_max). With this procedure, G_max and f_max are estimated again. This process is repeated twice.

Analysis of the parameter space exploration for cascades

The outputs of the five studies performed with cascades are analyzed in the following way. The criteria to identify signal termination is adapted from Section 1, being more flexible for cascades: levels 1 or 2 are said to have signal termination if the maximum of the temporal profile is higher than 0.1 and Q is higher than 1.2 (given that Q = P_max/P_ss, this condition means that there must be a decrease of at least 0.83% from the maximum to the steady-state to classify the output as having signal termination).

If P₁ and P₂ do not satisfy the criteria for signal termination, this set is labeled as monotonic for both cycles. If P₁ or P₂ (only one of them) satisfies the criteria for signal termination, then this set is labeled as signal termination in the first or second cycle, respectively. If, instead, both P₁ and P₂ satisfy the criteria, this set is labeled as signal termination in both cycles. Outputs, where P₁ and P₂ are not monotonic but do not reach the criteria for signal termination, are not considered for the analysis.

Frequency preference is evaluated in the cycle that exhibits signal termination, if only one of them does it, and in the second cycle if both of them have signal termination.

Procedure to obtain an adapted periodic input

The adapted sinusoidal input is defined as follows:

S_T = S_T,max (1 + 0.5sin (2πft))/1.5and the adapted train of square pulses is defined as follows:

S_T = S_T,max if mod(t,T) < d, S_T = 0 if not, where f is the frequency, mod means module, T is the period, and d is the length of the square. In this way, and given S_T,desired the sinusoidal input or the train of square pulses, we adapt S_T so that it is always positive. We do so by changing free substrate value S in this way: if S = S_Tdesired-P-C₁-C₂ < 0, then S = 0.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.