Properties of cell signaling pathways and gene expression systems operating far from steady-state

Ligand-receptor systems, covalent modification cycles, and transcriptional networks are basic units of signaling systems and their steady-state properties are well understood. However, the behavior of such systems before steady-state is poorly characterized. Here, we analyzed the properties of input-output curves for each of these systems as they approach steady-state. In ligand-receptor systems, the EC50 (concentration of the ligand that occupies 50% of the receptors) is higher before the system reaches steady-state. Based on this behavior, we have previously defined PRESS (for pre-equilibrium sensing and signaling), a general “systems level” mechanism cells may use to overcome input saturation. Originally, we showed that, given a step stimulation, PRESS operates when the kinetics of ligand-receptor binding are slower than the downstream signaling steps. Now, we show that, provided the input increases slowly, it is not essential for the ligand binding reaction itself to be slow. In addition, we demonstrate that covalent modification cycles and gene expression systems may also operate in PRESS mode. Thus, nearly all biochemical processes may operate in PRESS mode, suggesting that this mechanism may be ubiquitous in cell signaling systems.


Supplementary
Interval for simulation Interval from Huang and Ferrell [1] Interval from Bhalla and Iyengar [2] k 1 Table S1. For each of the parameters of the cycle, we indicate the interval considered for simulation and the intervals given in the paper by Huang and Ferrel and in that by Bhalla and Iyengar. This criterion for sampling intervals for a CM cycle was taken from [4].

Ligand-receptor system
The simplest binding model is a receptor that binds a ligand forming a complex , and is described by the following reaction: Using mass action kinetics and the conservation relation = + , the dynamics of the complex is described by the following differential equation: where R was obtained from the conservation relation, = − , to reduce the number of variables. We introduced the following definitions: ( ) = ( )/ is the amount of ligand relative to the dissociation constant for the binding-unbinding reaction ( = / ), = / is the amount of ligand-receptor complex relative to the total amount of receptors, and = / is the time expressed in units of a reference time = 1/ . x, y, and tn are dimensionless variables and transform equation (S2) into: We further assumed that the amount of ligand increases following an exponential function characterized by a maximum value and a characteristic time , resulting in = / and = / , which connects the time-scale associated with the ligand accumulation ( ) with that of the bindingunbinding process ( = 1/ ): We here consider two limit cases: fast stimulation with slow binding/unbinding and slow stimulation with fast binding/unbinding.
1a. Ligand-receptor system: fast ligand accumulation and slow binding/unbinding ( << 1) For << 1, equation S4 is reduced to ( ) ≅ , so this limit case resembles the system we described in [3] and in the Introduction of this paper, i.e, a receptor R interacting with a step-like temporal profile of ligand L. The solution of equation S2 is, then: The steady-state of equation S5, which is obtained by taking the limit → ∞, has an = 1 and 50 = 1. We are interested in obtaining the 50 and for versus , for any given time, which will result in two relationships: 50 ( ) and ( ).
To obtain these functions we need to solve the same equation with different values. Hence, we define a real number as the corresponding fraction: = 0.5 for the EC50, = 0.1 for the EC10 and = 0.9 for the EC90. Then, where is the value of the input that gives as the output (e.g. for = 0.5 , = 50 ). From equation S6 we can easily compute the time as a function of the remaining variables (i.e. the inverse function) Based on this last expression, EC50 versus time results in the plot we show in Fig. 1D in the main text. According to this plot, 50 decreases when time increases, reaching a value of 1 at binding equilibrium.
The Hill coefficient is defined as = 10(81)/ 10( 90 / 10 ) (S8) We obtained an approximated = 1.4 in the limit of low values of . This value is achieved by deriving ( 90 ) and ( 10 ) in the same way we derived ( 50 ), and neglecting terms as follows: These approximations are valid because at short times ( << 1), both 90 and 10 are much greater than one, similar to the behavior of 50 ( Figure 1D, main text). From the last two formulas we obtained the ratio 90 / 10 : Which results in 21.85 when is neglected, using equation (S8) this leads to = 1.42. This approximated result indicates that the sensitivity goes from about 1.4, in the limit of low values of to 1, in the limit of → ∞, and it means that before reaching equilibrium the sensitivity is higher than when equilibrium is reached.

1b. Ligand-receptor system: slow ligand accumulation with fast binding/unbinding
It is convenient, in this section, to re-define the dimensionless time so that it is expressed relative to the slowest time-scale in the system ( , ligand accumulation time). In this way: ∼ = / , resulting in where, as before, (~) = (~)/ is the normalized input and = / is the normalized output. = 1/( ) is a dimensionless parameter that is a ratio of the time characterizing the binding-unbinding reaction and the time characterizing ligand accumulation. We consider the limit of slow ligand accumulation with fast binding/unbinding, thus, ≪ 1. In this limit, the (binding) reaction is fast, equilibrating rapidly and remaining in near-equilibrium even as the variable slowly changes. Thus, we take the quasi-steadystate approximation / ∼ = 0, resulting in: The steady-state of the input-output curve is, as before, /(1 + ), a hyperbolic curve characterized by 50 = 1 and dynamic range 90 / 10 = 81. However, if we consider the time-dependent input-output curve, . ., versus for a given time ∼ , we find a time-dependent 50 : Interestingly, in this case the dynamic range is 81 (and thus = 1) for every value of ~.

1a. Mechanistic description and parameter values
A covalent modification cycle (CM) may be described by the following reactions: where is the kinase, is the phosphatase, S the unmodified substrate, S* the modified substrate, the complex between S and the enzyme that modifies the substrate (for example, if the enzyme were a kinase, by adding a phosphate group), and * the complex between S* and the enzyme that removes the modification (for example a phosphatase). ai are the association rates, di the dissociation rates, and ki the catalytic rates ( = 1 for the forward reaction and = 2 for the reverse reaction).
Applying the law of mass action, the kinetic equations governing the time evolution of such a system are: (S15.f) and, consequently, the conservation relations are We scanned parameter values randomly with log uniform distribution within the intervals defined on Table  S1 and following the approach described in the Methods section.

Kinetic rates for the examples included in Figure 3 B-I
For the Zero th order regime:

1b. Michaelis-Menten approximation
Using the Michaelis-Menten approximation and assuming that the total amount of substrate is ST=S+S* (in this sum we neglected the amount of substrate bound to the activating and deactivating enzymes), the system can be described by the equation: where x=(k1 Ea)/(k2 Ed) is the input to this system, expressed as the ratio of the maximum velocities of the activating and deactivating enzymatic reactions, with catalytic rate constants k1 and k2 ; y=S*/ST is the output, the fraction of active substrate; = , / and = , / are the Michaelis-Menten constants relative to the total amount of substrate; and time is expressed in units of a reference time tref=ST/(k2 Ed), so that tn=t/tref. .
As in the full mechanistic description, we assumed that the activating enzyme increases following an exponential function. This function is characterized by a maximum value Ea,max and a characteristic time Ea, resulting in: where = 1 , 2 and , = / .
Notably, simulation results in the four situations considered (first-or zero-order, fast or slow input) (Fig. S1) are in complete agreement with the results obtained with the full mechanistic description presented in the main text: 1) the input-output curve shifted from right to left over time; 2) when the stimulus increased slowly, the leftward shift in the input-output curve was correspondingly slow; 3) the shift is faster in zero-order, indicating that there is less time for PRESS when the enzymes are saturated; 4) while the nH decreased with time in the first-order regime, it increased in the zero-order regime, from about 1 towards its final high steady-state value (~25 in our simulation). Regarding the first-order regime, we noted in the main text that ( ) has a fast increase up to 1.42 and then decreases to 1. This initial increase up to 1.42, which does not appear in the Michaelis-Menten version, is probably due to the intermediary steps (complex formation): at short times, the overall rate of output production depends more on complex formation than catalysis itself.

Composing Shifters: A ligand-receptor activates a covalent modification cycle.
The full set of reactions for this model is composed of those for a ligand receptor (LR) system and those for a CM cycle, with an extra reaction involving the binding of the receptor and the substrate without the ligand (S19b). This last reaction does not lead to product formation. The reactions are as follows: For step stimulation, we described the system with the following differential equations: together with the conservation relations is the dissociation constant for the ligand-receptor reaction, is the total amount of receptor, is the unbinding rate for the ligand-receptor reaction, is the inactive receptor-substrate complex, is the ligand-receptor complex and the active enzyme for the CM cycle, the total amount of substrate, is the amount of free deactivating enzyme. * is the substrate in its active form and the substrate in its inactive form. The substrate binds the receptor in either of its two forms: free ( ) or bounded to the ligand ( ). Binding and unbinding rates for the substrate to R and RL are 1 and 1 , respectively. We assume that R binds L and S independently. Hence, the kinetic rates are identical. Finally, S* is produced and the enzyme ( ) released, with rate 1 . Similarly, * binds reversible to its deactivating enzyme with rates 2 and 2 , forming the intermediate complex * , and releasing and with rate 2 . Time is dimensionless, = / with = / ( 1 ) being a reference time. We define, as in the main text, the dimensionless Michaelis-Menten parameters = , / , and = , / , being , and , the dimensional values of the Michaelis-Menten constants.

Cascade of covalent modification cycles
The mechanistic model of a cascade of CM cycles is based on a three-tier mitogen activated protein kinase cascade, captured by the following reactions: To simplify the notation, we avoided MAP in each variable ( is named , and so on). Each of the three kinases has a dedicated phosphatase, termed after its substrate kinase with the suffix "P'ase". In this model, the reactants in abundance (ATP, water, and so) are assumed to be constant, so they are included in the rate constants.
The 10 reactions described above give rise to 18 differential equations: The reaction rates are only given as and in the literature. Thus, we estimated (enzyme binds substrate), (enzyme releases substrate without modifying it) and (enzyme modifies and releases the substrate) according to the Michaelis-Menten formula and the additional assumption, that and have a constant ratio = / . Then the rate constants can be calculated by: = , , = (1 + ) , = , . As proposed by [5], we choose = 4. We scanned 100 values of the input 1 from 10 −3 for 1 in logarithmic scale for different times, the output was the concentration of the doubly phosphorylated MAPK,[ ]. We computed Input-Output curves using the ode23s function from Matlab to integrate.

Time window and maximal EC50 correlation
Analyzing the values of TW and ME for each parameter set, we found that they have a high correlation in the log-log space, ( 10( ), 10( )) = = 0.8562.
In the main text we showed the results for the Maximal EC50 for all combinations of parameter pairs. Due to the high correlation just described, the heatmaps for the Time Window are qualitatively similar, and leads to similar conclusions ( Figure S2).