The Goldbeter-Koshland switch in the first-order region and its response to dynamic disorder.

In their classical work (Proc. Natl. Acad. Sci. USA, 1981, 78:6840-6844), Goldbeter and Koshland mathematically analyzed a reversible covalent modification system which is highly sensitive to the concentration of effectors. Its signal-response curve appears sigmoidal, constituting a biochemical switch. However, the switch behavior only emerges in the 'zero-order region', i.e. when the signal molecule concentration is much lower than that of the substrate it modifies. In this work we showed that the switching behavior can also occur under comparable concentrations of signals and substrates, provided that the signal molecules catalyze the modification reaction in cooperation. We also studied the effect of dynamic disorders on the proposed biochemical switch, in which the enzymatic reaction rates, instead of constant, appear as stochastic functions of time. We showed that the system is robust to dynamic disorder at bulk concentration. But if the dynamic disorder is quasi-static, large fluctuations of the switch response behavior may be observed at low concentrations. Such fluctuation is relevant to many biological functions. It can be reduced by either increasing the conformation interconversion rate of the protein, or correlating the enzymatic reaction rates in the network.

Introduction

A biological system usually functions by regulating protein activities through protein interaction networks (PINs), which are built of interconnected modules. Some typical modules were summarized in a recent review by Tyson et. al.[1]. One of the most common modules in the PINs is the covalent modification system, which typically consists of a phosphorylated-dephosphorylated couple: [2], [3], [4]. Here R is the protein being covalently modified; its phosphorylated form, R, amounts to the response of the system. The kinase, S, enters as the external signal. A denotes the phosphatase, which restores the protein to its ‘non-response’ form. This system shows “zero-order ultrasensitivity”: sharp transitions occur in the signal-response curve when the modification enzymes, S and A, are saturated by the substrate, R and R, i.e. [S] and [A] much smaller than [R] = [R]+[R]. Eq. (1) gives the steady-state response of the system [4],where K 1 and K 2 are the Michaelis-Menten constants. According to Eq. (1), [R] appears as a sigmoidal, Goldbeter-Koshland function of the kinase concentration, [S]. The system thus behaves like a switch in response to external signals, constituting an important module in the PIN [5], [6].

However, the assumption of saturated enzyme reaction does not always hold in real PINs. When [S] and [A] become comparable to [R], the system transits from the zero-order into the first-order region, and the ultrasensitive switching behavior disappears in the simple two-component system. In the first part of this paper, we will show that the switching behavior in the first-order region can be restored by an additional cooperative mechanism of the phosphorylation reaction.

In the second part, we will discuss the effect of dynamic disorder on the biological switch. Dynamic disorder refers to the phenomenon that the ‘rate constant’ of a reaction appears as a stochastic function of time [7], [8].This phenomenon has attracted extensive experimental and theoretical studies since the pioneering work on ligand binding to myoglobin by Frauenfelder and coworkers [9]. Recently, single-molecule experiments directly showed the dynamic disorder in enzymatic reactions [10], [11], [12]. We refer the readers to some review articles and references therein for more information [13], [14], [15], [16].

The dynamic disorder results from the relatively slow fluctuations of the protein conformation, either in the enzyme or the substrate. Conformation fluctuations result in fluctuations of the enzymatic reaction rate. In the traditional chemical reaction theories, it is assumed that the reaction occurs much more slowly than the conformational fluctuations; so the reaction rate observed on the slow reaction time scale appears as the ensemble average of the reaction rates of each conformation. But biochemistry and biophysics studies showed that the conformational fluctuations of proteins can be much slower than previously assumed, because of the rugged energy landscape of protein conformations [10], [11], [13], [17], [18], [19]. The fluctuation time scale is sometimes comparable or even slower than that of the reaction itself. In this regime, the observed reaction rate becomes stochastic in time, reflecting the fluctuation of the enzyme conformations.

In this work, we will show that in the switch module studied in the first part, the dynamic disorder induced by the conformational fluctuations of the substrate protein mainly affects the variance of the system response, but not the ensemble average response. We also investigated two ways to reduce the noise originated by the dynamic disorder.

Results

The sigmoidal switch outside the zero-order region

The covalent modification system only achieves zero-order ultrasensitivity when the substrate proteins are much more than the enzymes. Without this assumption, the intermediate products should not be reduced from the reaction pathway. The full pathway is presented in case a, Figure 1. The signal-response (SR) relationship of the system is given by the steady state solution of the governing equations (Eq. (2) and (3)):with concentration constraintsFor mathematical simplicity, we assumed that the phosphatase is in such great excess that its concentration, [A], remains approximately constant throughout time. So we absorbed it in k 2. Similar treatment was made throughout the paper. Numerical calculations confirmed that relaxing this assumption does not qualitatively change the conclusion of this paper (not shown).

Figure 1

Candidate schemes for the sigmoidal switch.

Case a (upper panel within the shaded box): the full pathway of the zero-order ultrasensitive switch studied by Goldbeter and Koshland. The protein R is phosphorylated and dephosphorylated by enzyme S and A respectively. The phosphorylated form R represents the response and the enzyme S the signal. Case a2 (the entire upper panel): R and S form a loosely bound complex, then a tightly bound one preceding the phosphorylation step. Case b (lower panel, the shaded area is the same as case a): RS can bind another S molecule, which accelerate the phosphorylation reaction k 1.

Experimentally, one can control as the external signal either [S], the total substrate concentration, or [S], the free substrate concentration (the d/dt[S] equation is not needed in this case). Unfortunately the above system does not produce a desired sigmoidal SR curve with either signal form. A sigmoidal SR curve must have zero second derivative of the response to the signal at the inflection point. But Eq. (2) and (3) give identically negative second derivatives of [R] to both [S] and [S] (relation (4), derived in Appendix S1, case a), which indicates the lack of sigmoidal behavior. The numerical result of this case is shown in Figure 2a.Adding a tight-binding step to the phosphorylation reaction does not produce a sigmoidal curve, either. In case a2 (Figure 1), for example, R and S first form a weakly bound compound, RS. RS then convert to the tightly bound form, R*S. Eventually R*S proceeds to phosphorylation. In this system, the second derivative of the response to the signal is also mono-signed (see Appendix S1, case a2). Therefore no sigmoidal behavior emerges.

Figure 2

SR curves of the candidate schemes.

(a) Example SR curves obtained with the schemes shown in Figure 1. The G–K function is the result Goldbeter and Koshland obtained for the zero-order region. Parameters of other cases are obtained by fitting the steady state solutions to the G–K function (see Methods), and are listed in Table 1. (b) Concentrations of various protein forms as a function of the total signal concentration in case b, the only scheme showing the switch behavior.

Table 1

Model parameters.

	k 1	k 1f	k 1r	k 1f′	k 1r′		k 2	k 2f	k 2r	[R]t
Case a	5	0.06	20				1	0.06	20	1
Case b	0.006	400	12	40	670	10	1	4	79	1

The above analysis suggests that nonlinear terms of [S] are required to generate the sigmoidal behavior. We examined one such scheme (case b, Figure 1), inspired by the work of Sabouri-Ghomi et al. [20]. In this case, we modified the model such that binding an additional S molecule to the intermediate compound, RS, greatly facilitates the phosphorylation (k 1′≫k 1 in Eq. (5)). The governing equations becomewith concentration constraintsThe additional mechanism brings about a nonlinear term of [S] (see Eq. (5)) and leads to the desired sigmoidal response (Figure 2a). Figure 2b shows the concentrations of the various protein forms as functions of the signal [S]. At low kinase concentration, most protein substrates are sequestered in the single-kinase intermediate, RS; in this regime, [RS] increases steadily with [S]. Beyond a critical kinase concentration, though, the double-kinase intermediate, SRS, starts to form and stream through the phosphorylation reaction. This turnover results in the corresponding sudden drop in [RS] and abrupt rise in the final signal, [R]. Our numerical analysis also revealed bistability in this case, which was confirmed by the analysis with the Chemical Reaction Network Toolbox [21]. This scheme of bistability is additional to that reviewed by Kholodenko [6]. Further analysis is necessary for the bistability.

For the same circuit, we also carried out the calculation when the control signal is the free kinase concentration, [S] (reducing the d/dt[S] equation). Our numerical analysis showed that the SR curve, i.e. the [R]−[S] curve, is only mildly sigmoidal over a wide parameter range.

Effects of dynamic disorder

In this part, we will investigate the effect of dynamic disorder on the sigmoidal switch, in particular, on the circuit presented in case b from the previous section. We adopted a multistate model studied by Kou et al. [22], in which the substrate protein assumes N different conformations, R, …, R. All the other forms of the protein possesses N corresponding conformations, for example, RS,…, RS for RS. Only matching conformations of the reactant and the product are admissible in chemical transitions, e.g. R→AR, but not R→AR. Some or all of the chemical reaction rates vary with different conformations. As the proteins randomly change their conformations, the average rates of the reactions undergo temporal fluctuations, or dynamic disorder. Since loops of reversible reactions exist in the system, e.g. R↔R↔RS↔RS↔R, the reaction rates on the loops have to satisfy the constraint of detailed balance. The products of rates in the two directions of the loop have to be equal. To save the trouble of maintaining this constraint, we only considered the dynamic disorder in the irreversible reactions, like SRS→R, i.e. only k 1,k 1′,k 2 are allowed to fluctuate along the conformational coordinate. Interconversion occurs between each pair of conformers of the same protein form. For simplicity, a uniform rate was used for all the conformation interconversion. The governing equations of the system consist of Eq. (5) repeated over the N conformations, as well as the equations for conformation interconversion.

This model can be regarded as a discrete representation of the continuous, coupled diffusion-reaction model widely used in studies of protein motors and other macromolecules [23], [24], [25], [26], [27], [28], [29], [30]. In the continuous model, the molecule of interest assumes several chemical states. At a given chemical state, the molecule diffuses along a conformational coordinate. Transitions between the chemical states happen vertically, i.e. without simultaneous displacement along the conformational coordinate. This is because the chemical transition is generally characterized as a barrier-crossing process: the waiting time for a transition (determined by the rate constants) may be long, but the actual transition happens on a very short time scale. This time scale usually allows no resolvable displacements along other degrees of freedom. In our model, the different complexes (RS, SRS, etc.) are analogous to the chemical states in the continuous model, and the conformational coordinate is discretized, accounting for the energy minima on the rugged energy landscape.

We first studied the dynamic disorder of k 1′, since the facilitated phosphorylation, SRS→R, is the key step that restores the desired switching behavior in case b (Figure 1b). We assumed that k 1′ obeys the gamma distribution over the conformations (mathematical expression given in Methods). Figure 3 compares the resulting ensemble SR curve, vs, with the SR curve of each conformation, [R]vs[S]. Here and are the ensemble concentration of the signal and the response, i.e. , . And [S] is the scaled signal concentration for the i-th conformation, [S] = [S]+N([RS]+2[SRS]). In the limit of zero conformation interconversion rate, the SR curve of the i-th conformation corresponds to the scenario as if only this conformation exists.

Figure 3

SR curves of case b in the presence of dynamic disorder in k 1′.

A set of k 1′ of 10 different conformations are computed from a gamma distribution p(k) = [1/(bΓ(a))k −1exp(−k/b) (see Methods). For the results shown here, a = 4, b = 10/a, so that the mean rate constant is 10, the k 1′ value used in the first section (Table 1). (a) The ensemble averaged response curves with the conformation k 1′ interconversion rate r int = 10−5 (solid line) and r int = 10−1 (dashed line). (b) Responses by individual conformers with r int = 10−5 (dotted lines). The ensemble averaged response is also shown in comparison as the solid line. (c) Same as b with r int = 10−1.

Figure 3a shows that the ensemble SR curve is not significantly affected by the conformation interconversion rate, r int. Therefore, the effect of dynamic disorder on the system is barely noticeable if the response is measured at bulk. But the variance of the response strongly depends on r int. When the conformation interconversion is much slower than the chemical reactions, the reflection point of the SR curve shifts up to 20% of the original value (Figure 3b). This variance is dramatically diminished when r int increases to two orders of magnitude below the mean value of the disordered reaction rate (Figure 3c). Figure 4 summarizes the variances resulting from several trial cases. Dynamically disordered k 2 (case a) generates comparable variances as the disordered k 1′ (case b) does, while disordered k 1 (case c) generates much smaller variances. This is because k 2 and k 1′ are both associated with the main reaction pathway in the system, but k 1 is not (cf. Figure 1b). Case d of Figure 4 shows that perfect correlation between the disordered rates can also diminish the variance. In this case, k 1, k 1′ and k 2 share exactly the same distribution over the conformations; the variance of the response reduces and becomes invariant to the interconversion rate, r int. This happens because change of the kinase activity is perfectly compensated by the change of the phosphatase activity. However, this phenomenon is not robust at all. As shown in case e of Figure 4, broadening the distribution of k 1′ immediately restores the large variance.

Figure 4

Variances of the SR curves change with the conformation interconversion rate, r int, with different sets of disordered parameters.

Upper: r int vs. relative variances of the critical signal level. The critical signal level is defined as [S] at half the plateau responses (i.e. [R] = [R]([S]→∞)/2). Bottom: r int vs. relative variances of the plateau response (i.e. [R]([S]→∞). Different sets of disordered parameters denoted in the legend: a) k 2 is computed from the gamma distribution with a = 4, b 2 = 1/a; b) k 1′ is computed from the gamma distribution with a = 4, b 1′ = 10/a; c) k 1 is computed from the gamma distribution with a = 4, b 1 = 0.008/a; d) all the enzymatic reaction rates, k 1, k 1′ and k 2, are computed from the gamma distribution with (a,b 1;a,b 1′;a,b 2); e) and k 2, come from the gamma distribution with (a,b 1;a,b 2), but k 1′ comes from the gamma distribution with (a/2,2b 1′). The parameters of the gamma distributions in case a, b, c are chosen so that the mean rates equal the ones used in absence of the dynamic disorder. Note that the parameter a, which determines the width of the distribution, were chosen the same except in the last case for k 1′.

Discussion

This work has two focuses. First we discussed the detailed mechanisms of the switch module in the first-order region. It is well known that the covalently-modification system studied by Goldbeter and Koshland cannot produce the switching behavior in this region. To recover this behavior, we grafted onto the Goldbeter-Koshland system the cooperative binding mechanism, which produces sigmoidal responses by itself. When the reaction is facilitated by the binding of two enzyme molecules, the system regains the remarkable Hill factor in the first-order region. While experimental studies on the proposed mechanism is lacking, it is common for enzymes to work in the form of dimers during signal transduction. Further experimental studies are necessary to test the theoretical result for the basic switch module here.

Next we discussed the effect of dynamic disorder on the switch module. The ensemble averaged behavior is nearly insensitive to dynamic disorder. But when the number of protein molecules is small, fluctuations of the module response will significantly affect the functioning of the switch. The critical signal level for the transition of response can shift up to 20% when the reaction rates fluctuate slowly. This study raises two important questions. First, has the nature evolved to reduce the effect of dynamic disorders in PINs, especially some vital ones that require high robustness? Our study suggests two possible mechanisms to reduce the noise: increasing protein conformation interconversion rates, or correlating the distributions of the chemical reaction rates in the PIN. Second, can the dynamic disorder in one module of the PIN actually reduce the overall noise of the PIN? Noise effect and reduction in a biological network is an actively studied field [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], [43], [44]. The dynamic disorder, as one source of noise, can have broad time scales, and interplay with other noises. Theses noises may offset each other in the overall behavior. For example, in this work only the substrate protein has different conformations. But in a real PIN with cascades of such phosphorylation-dephosphorylation cycles, the substrate of one cycle serves as the kinase or phosphatase of the next cycle. Novel behaviors will probably emerge with the complete set of dynamic disorders built in. This is a future direction to follow.

In this work we focused on the phosphorylation-dephosphorylation cycle. Same conclusions apply to other signal transduction modules with similar kinetic structures. One example is the GTPase cycle [45], [46], [47]. A membrane-bound GTPase switches between a GTP-bound form, and a GDP-bound form, which assume different enzyme activity. Such a protein can detect GTP concentration as the signal, and respond with effects on the downstream biochemical pathways. The dynamics of such GTPase cycle is mathematically isomorphic to that of the phosphorylation-dephosphorylation cycle [47]. Thus, it should demonstrate the same switching behavior and effect of the dynamic disorder.

Methods

In the first part, we used both the downhill simplex method and the simulated annealing method to search the parameter space of rate constants for the minimum of the following function [48]:Here [R(S)] is the calculated R concentration as a function of the signal strength, and T(S) is the desired sigmoidal SR curve. In all the calculations T(S) assumes the ‘Goldbeter-Koshland’ function form,with u = 2S, v = 1, J = K = 0.05.

In the second part, the dynamically disordered enzymatic reaction rates were computed from a gamma distribution with mean ab and variance ab (Eq. (9)). a and b are given in the corresponding figure captions.The above continuous distribution is discretized in the following way. First the parameter coordinate was divided into N bins with equal accumulated probability, i.e.where ξ's are the boundaries of the bins, with ξ0 = −∞, ξ = ∞. Then the N discretized rate constants were chosen as,

(0.11 MB DOC)

Click here for additional data file.