Originally published as Biophys J. BioFAST on February 24, 2006.
doi:10.1529/biophysj.105.060491
Biophysical Journal 90:3749-3761 (2006)
© 2006 The Biophysical Society
Regulated Degradation Is a Mechanism for Suppressing Stochastic Fluctuations in Gene Regulatory Networks
Hana El-Samad and
Mustafa Khammash
Mechanical and Environmental Engineering, University of California at Santa Barbara, Santa Barbara, California
Correspondence: Address reprint requests to M. Khammash, Tel.: 805-893-4967; E-mail: khammash{at}engineering.ucsb.edu.
 |
ABSTRACT
|
|---|
Cellular events that execute life programs are ordered and reproducible, despite the noise and randomness underlying their biochemical reactions. The identification of the processes that ensure this robust operation is essential not only to uncover the salient design principles in organisms, but also to forward-engineer reliable genetic networks for biotechnological and therapeutic purposes. The use of feedback for noise reduction has been suggested as a recurring motif in genetic networks. In this work, we show that regulated degradation of proteins implements a negative feedback loop that enhances robustness against stochastic fluctuations and cellular noise. The analysis is carried out in the context of the bacterial heat shock response where the tight control of the amount of heat shock proteins is achieved through an intricate architecture of feedback loops involving the 
32-factor. The 32 regulates the transcription of heat shock proteins under normal and stress conditions. An essential feature of the heat shock response is a feedback loop regulating the degradation of 32. We investigate the noise-rejection properties of this loop, therefore illustrating our point in a biologically plausible example.
|
INTRODUCTION
|
|---|
Living systems are inherently noisy. In addition to their exposure to environmental noise, they are subjected to biochemical noise that leads to stochastic fluctuations in their molecular species. The magnitude and nature of the fluctuations is thought to be dependent on the network structure, the concentration of molecules that populate this structure, and the reaction rates of the underlying biochemical reactions. At the same time, these systems are known to function reliably and even thrive in the presence of noise. It has been suggested that this robust operation is the result of intricate noise attenuation and exploitation mechanisms, such as feedback regulatory loops that either control the concentration of important effector molecules (1
,2) through regulated synthesis or their activity and availability through multimerization and sequestration (3–5). Coupled with other mechanisms that implement high Hill coefficients and time lags, feedback loops have been shown to yield dynamics that exhibit sophisticated behavior, including phenomena of noise suppression by noise (6,7).
One of the earliest theoretical investigations of the role of feedback in noise attenuation was given in Thattai and van Oudenaarden (1). This study (1) involved a simple model of prokaryotic transcription. Two scenarios were considered: an unregulated gene and a gene where the protein regulates its own production in a negative feedback fashion. For the linear model used to describe the system, the statistical moments of the various species are solvable analytically from the chemical master equation. It was found that the autoregulatory system has a Fano factor (the variance of the protein steady-state distribution normalized by the mean value of the distribution) which is closer to unity, indicating smaller deviation of the protein distribution from Poissonian statistics in the case of feedback. Experimental evidence of this role of protein synthesis was partly addressed in Becskei and Serrano (8).
Negative feedback is also frequently implemented through the control of protein degradation. In the heat shock system of Escherichia coli, for example, the stability of the heat shock factor 32 is tightly controlled by an autoregulatory loop whereas its synthesis is directly influenced by the ambient temperature and is only feedforward-regulated. In addition to stress responses, both regulated protein and mRNA degradation play a crucial role in a large number of cellular functions such as cell cycle regulation, growth, division, and differentiation (9–13). In fact, protein proteolysis is a necessity for maintaining homeostasis while cellular components are continually rebuilt in response to external inputs or during development. Furthermore, protein degradation provides a mean to terminate the life of regulatory proteins, such as cyclins, transcription factors, and components of signal transduction pathways in a timely manner (14), providing the prospects for intricate regulation. For instance, the RNA binding protein and global regulator Hfq whose disruption causes decreased growth rates and osmosensitivity is known to undergo degradation-based translational control (12). On the protein level, the controlled ubiquitin-mediated degradation of tumor suppressors and proto-oncogenes forms the basis of healthy cell growth and proliferation. An extensively studied example is that of the tumor suppressor p53 that transcriptionally activates Mdm2. Mdm2 in turn targets p53 for degradation (15). A similar logic is employed for the controlled proteolysis of many transcription factors such as the nuclear factor 
B (NF–B), a ubiquitous inducible transcription factor involved in the central immune system, and the inflammatory, stress, and developmental processes (16).
In this article, we present a stochastic study of the properties of regulated degradation in the heat shock response system. We first demonstrate the noise attenuation advantage conferred by this strategy in a simple network where the end product is involved in the regulation of the stability of its precursor. Since a negative feedback loop can also be implemented in the same example when the end product regulates the synthesis of its precursor, we offer insight into some similarities and differences between the two cases. We then present a more biologically relevant example involving the heat shock response where regulated proteolysis plays a crucial role. Indeed, regulation of the degradation of the -factor, 32, is an intriguing feedback mechanism by which the heat shock response in the bacterium E. coli is partly controlled. This degradation effect influences various performance aspects of the heat shock response. For example, this loop enhances the transient response to a sudden increase in the temperature in the cell (17). We show here that this loop has an additional crucial benefit. It alleviates the harmful effects of stochastic fluctuations, therefore producing a more reliable chaperone signal to tackle the delicate protein folding process. Our analysis is based on the stochastic simulation algorithm (SSA), which predicts a narrower distribution for the number of cellular chaperones as a result of degradation feedback. In the limit of small noise, and in terms of the linear noise approximation (LNA), this result is attributed to an increase in the feedback gain of the regulated system as compared to schemes where degradation is constitutive.
|
NUMERICAL AND ANALYTICAL METHODS FOR STOCHASTIC CHEMICAL KINETICS: THE STOCHASTIC SIMULATION ALGORITHM AND THE LINEAR NOISE APPROXIMATION
|
|---|
In principle, the exact stochastic behavior of the populations of different species involved in genetic models is described by the chemical master equation (CME), a differential equation governing the evolution of their joint probability density function. The CME describes the time evolution of the probability of having a certain number or concentration of molecules, as opposed to deterministic rate equation descriptions of the absolute concentration of these molecules (18,19). In the CME, reaction rates are transformed into probability transition rates which can be determined based on physical considerations. The CME can be derived based on the Markov property of chemical reactions. Using this Markov property, one can write the Chapman-Kolmogorov equation, an identity that must be obeyed by the transition probability of any Markov process. Using stationarity and taking the limit for infinitesimally vanishing time intervals, one obtains the Master equation, as the differential form of the Chapman-Kolmogorov equation (18,20). The CME can be written fairly easily from the rate constants and stoichiometries of the reactions in a chemical system. However, more often than not, it can neither be solved analytically nor numerically. Therefore, one has to resort to a Monte Carlo type of numerical simulation, such as the Gillespie stochastic simulation algorithm (SSA) (21). In general, one needs to break the biochemical interactions in a system into elementary reaction steps (birth and death processes) characterized by their probability of occurrence (otherwise called propensity functions). Based on these propensity functions and at any given time instant, one can compute the type of reaction that will occur next in the system as well as the time at which such a reaction would take place. Iteration of this procedure until a final time is reached provides numerical realizations possessing the same statistical properties described by the master equation.
A less computationally demanding, yet approximate, approach is the simplification of the master equation in a linear noise approximation (LNA). Roughly speaking, the LNA involves the expansion of the master equation in a Taylor series near macroscopic system trajectories or stationary points. Terms of first order in the expansion are then identified with the macroscopic rate equations, whereas terms of second order describe the approximate noise acting on the system. More specifically, approximate explicit expressions for the stationary covariance matrix Cs containing the variances and covariances of the fluctuations in the components of the system is obtained as the solution of an algebraic Lyapunov equation,
 | (1) |
As is the Jacobian matrix of the deterministic reaction rate equations evaluated at the steady state, i.e., 
where S is the stoichiometry matrix, W is the propensity vector, and 
s is the deterministic steady-state solution. The value Ds is the diffusion matrix that embodies the contribution of the chemical reactions to the overall fluctuations. In particular, Ds takes the form Ds = S[diagW(s)]ST, where S and W are as defined above and [diagW(s)] is a square matrix having the elements of W(s) on the diagonal, and zero elsewhere. The LNA was known in the early literature as the van Kampen's system size expansion (18) and was recently rederived for multivariate CMEs in Elf and Ehrenberg (22) and Tomioka et al. (23).
In our analysis, we shall exclusively use the results of the SSA to make quantitative and qualitative conclusions about the various systems under consideration. In the limit of small noise, we complement these results with the more analytical, but inevitably approximate results of the LNA. The LNA formulation can potentially single out, at least qualitatively, the effects of the various structural mechanisms of a system (such as feedback loops) and their role in noise attenuation or amplification. The validity of these predictions is, however, limited by the incapability of the LNA to capture some inherently nonlinear behaviors. Whenever this is the case (as will be shown through an example in Regulated Synthesis Versus Regulated Degradation, in the section below), one needs to resort back to the results of the SSA.
|
REGULATION OF PROTEIN PROTEOLYSIS IS A MECHANISM FOR NOISE ATTENUATION: A SIMPLE EXAMPLE
|
|---|
In this section, we present a simple example where the essence of the regulated proteolysis mechanism can be extracted and analyzed. In this system, a protein X1 is produced from a constant pool of substrate. The protein X1 promotes the production of another protein X2, which is then degraded with first-order kinetics at a rate 
2. We consider two scenarios for the degradation of protein X1. In the first scenario, X1 is degraded at a constitutive rate 0. We will refer to this case as an open-loop scheme. In the second scenario, protein X2 regulates the degradation of X1 in a closed-loop fashion. The degradation of X1 should therefore appear as a function of x2 (we use lower case letter x to denote the number of molecules of X), f1(x2) multiplied by x1. For simplicity, we choose a linear function f of x2, i.e., f1(x2) = rx2. The two regulation scenarios are illustrated in Fig. 1, a and b. Every event involving the production or disappearance of the species x1 and x2 is a molecular reaction which can occur with a propensity that can depend on the kinetic rate constants and abundance of x1 and x2. The list of these chemical reactions, their propensities, and the resulting stoichiometric changes is given in Table 1.
View larger version (9K):
[in this window]
[in a new window]
|
FIGURE 1 Simple circuit that illustrates the use of negative feedback for noise attenuation. (a) Open-loop system. (b) Closed-loop system with regulated proteolysis. (c) Closed-loop system with regulated synthesis.
|
|
Exact numerical analysis using the SSA
Using the SSA and the elementary reactions in the Table 1, we simulated the two systems implementing constitutive and regulated degradation of X1. The parameters were randomly chosen except for the requirement that 
where we have used the superscript s to denote the steady-state number of X2. The reason for this choice is the necessity to implement the same mean values of X1 and X2 for direct comparison between the two schemes. Further, it corresponds to the physically relevant situation where two regulators maintain the same concentration in the cell but differ in the details of how it is achieved. In the linear degradation feedback case, this would correspond to setting 
The results of 250,000 steady-state samples are compiled as histograms for the numbers of the end-product of the network, protein X2. These histograms are shown in Fig. 2 a for the open-loop case and in Fig. 2 b for the closed-loop case. It is noticeable that the regulated proteolysis case is characterized by a distribution which is more peaked, and less spread around its mean. In the constitutive proteolysis case, on the other hand, x2 is more likely to take rather large excursions around its mean. Next, using the LNA, we analytically investigate the roots of this behavior for a general set of parameters.
Analysis using the LNA
To use the LNA, we form the stoichiometry matrix S and the propensity vector W(x1, x2) of the system as
and
where f(x2) = rx2 for the regulated proteolysis case and f(x2) = 0 for the constitutive proteolysis case. The steady-state Jacobian for the open-loop system in Fig. 1 a, denoted here by Acs, is given by
while the stationary Jacobian of the closed-loop system in Fig. 1 b, denoted by As, is given by
Now, we ensure internal consistency between the two models by requiring that all the variables assume the same steady-state values in the open-loop and closed-loop schemes. As in the previous section, all other parameters being equal, we ensure that x1 is degraded at the same rate, and hence assumes the same steady state in both models by setting 0 = 
in the open-loop system. Hence, 
s = cs, s = cs, and 
s = cs. The two systems possess diagonal diffusion matrices. These diffusion matrices are identical because the steady states of the systems variables are equal. We denote the nonzero terms of these diagonal matrices by d11 and d22. Using all the quantities above, we can now solve the algebraic Lyapunov Eq. 1 and extract the variance in the level of protein X2 for the two cases. The variance for the open-loop (constitutive degradation) system is given by
while that for the regulated degradation case is given by
The difference can be computed as
Notice that the difference between 
is always a positive and monotonically increasing function of 
S, the feedback gain. Therefore, the regulated degradation indeed possesses a lower variance than the constitutive degradation. This result confirms the intuition provided by the exact SSA output, and points to the extra feedback gain provided by the dynamic regulation of stability as the structural component responsible for the enhanced noise attenuation property.
Regulated synthesis versus regulated degradation
For the transcription module considered in the previous example, negative feedback can also be implemented by a negative feedback from X2 onto the synthesis rate of X1, as shown in Fig. 1 c. In this case, instead of a constant ß, the synthesis rate of X1 is given by f2(x2) where, for best comparison with regulated degradation, we choose 
The degradation of X1 is again constitutive and its propensity function is given by sx1. With this choice of f2, the same steady-state and average lifetime for x1 and x2 can be simultaneously achieved by setting 
and 
The parameters dictating the behavior of X2 are still given by the second column of Table 1.
The linear noise approximation applied to the transcription module with this synthesis feedback in place generates Jacobian and diffusion matrices identical to those of regulated degradation. This in turn implies that these two mechanistically different schemes of regulation generate the same variances for X1 and X2 in the case where fluctuations around the equilibrium are small enough to warrant the validity of the LNA. Underlying this conclusion is also the assumption that the mean of the distribution generated by either regulated synthesis or degradation is adequately represented by the deterministic steady state, and that the resulting distributions are unimodal. This is certainly the case for some regions of the parameter space. For example, for ß = 100, r = 1, k = 1, 2 = 1, the deterministic steady state generated by the differential equation description of the system assumes exactly the value 10. In this case, the open-loop and both regulated degradation and synthesis systems generate pdfs for X2 whose mean values are close to 10 and very close to each other. Furthermore, the noise attenuation shown by regulated synthesis and degradation are very similar (see Fig. 3 a), except that regulated synthesis occasionally assumes tail values that are not reached by regulated degradation. However, this cannot be generalized to situations where the inherently distinct nonlinear behavior of the two schemes dictates not only the magnitude of the fluctuations, but also the response of the system to them beyond what is captured by the LNA. When ß = 0.01, r = 0.0002, and k = 0.2, 2 = 0.01, the inherent nonlinear properties of the two schemes of regulation translate into differentiating stochastic properties (Fig. 3 b), the most striking being a pdf for X2 showing a bimodal shape for regulated degradation. The pdf has a small peak at X2 = 0, which is by no means a numerical artifact. In fact, one of the nonlinear characteristics of these models of regulated degradation and synthesis are trajectories that decrease before increasing when X1 and X2 simultaneously drop below their equilibria. However, in the phase that precedes the increase, regulated degradation exhibits a larger overshoot than regulated synthesis in the decreasing X2 direction. During this typical overshoot, the trajectory of X2 given by regulated degradation can hit the value zero while that given by regulated synthesis, for the same parameters and perturbation, cannot. Furthermore, when this happens, and in addition X1 is very small, the nature of the regulated degradation scheme is such that X2 has to wait for X1, which accumulates at a maximal rate ß, before it can increase again. This offers an explanation for the small peak at zero. In contrast, a small value of X2 makes the synthesis rate of X1 increase sharply since it is dictated by 
generating more opportunities for X to correct fast for the deficiency in its value. In fact, X2 does not hit zero and hits X2 = 1 only twice in 820,000 samples.
In addition to a small peak characterizing regulated degradation, the pdf is also characterized by a larger peak close to the deterministic steady state (the deterministic steady state is 
and the mean is 
and the higher peak occurs at 
). Interestingly, the regulated synthesis system simulated using either 
generates pdfs that are minimally different from each other, whose mean is much larger than 
and whose upper tail is larger than that of regulated degradation. In parallel, the pdf generated by regulated degradation has a substantially fatter lower tail at low molecule values.
Note that in this case, the open-loop system still shows a pdf that is less peaked and more spread than both regulated degradation and synthesis, again confirming our intuition about the role of feedback in noise attenuation and motivating the investigation of the role of regulated degradation in the natural setting of the heat shock response.
The heat shock response in E. coli: deterministic and stochastic models
High temperatures cause cell constituents such as proteins to undergo structural changes: proteins unfold from their normal shapes, resulting in malfunctioning and eventually death of the cell. Cells have evolved gene regulatory mechanisms to counter the effects of heat shock by expressing specific genes that encode heat shock proteins whose role is to help the cell survive the consequence of the shock. Many heat shock proteins serve as molecular chaperones that assist in the refolding of denatured proteins; others are proteases that degrade and remove the denatured proteins. In E. coli, the heat shock (HS) response is implemented through an intricate architecture of feedback loops centered around the heat-shock factor (a molecule referred to as 32) that regulates the transcription of the HS proteins under normal and stress conditions. At physiological temperatures (30–37°C), there is very little 32 present and hence little transcription of the HS genes. When bacteria are exposed to high temperatures, 32 first rapidly accumulates, allowing increased transcription of the HS genes and then declines to a new steady-state level characteristic of the new growth temperature. An elegant mechanism that senses temperature and immediately reacts to its effect is implemented as follows in the bacterial heat shock response. At low temperatures, the translation start site of 32 is occluded by basepairing with other regions of the 32 mRNA. Upon temperature upshift, this basepairing is destabilized, resulting in a melting of the secondary structure of 32, which enhances ribosome entry, therefore increasing the translation efficiency. Indeed the translation rate of the mRNA encoding 32 increases immediately upon temperature increase (24). Hence, a sudden increase in temperature, sensed through this mechanism, results in a burst of 32 and a corresponding increase in the number of heat shock proteins. This mechanism implements a control scheme very similar to a feedforward control loop typically utilized in many engineering systems. With this mechanism in place, the production of heat-shock proteins becomes temperature-dependent. In addition, two feedback mechanisms are also used in the regulation of the activity and degradation of 32 as follows. The chaperone DnaK and its co-chaperone DnaJ perform the function of protein folding. At the same time, they are capable of binding to 32, limiting its capability to bind to RNA polymerase. Increased temperature produces an increase in the cellular levels of unfolded proteins that then titrate DnaK/J away from 32, allowing it to bind to RNA polymerase and resulting in increased transcription of DnaK/J and other chaperones. The accumulation of high levels of HS proteins leads to the efficient refolding of the denatured proteins, thereby decreasing the pool of unfolded protein, freeing up DnaK/J to sequester this protein from RNA polymerase. This implements what we shall refer to as a sequestration feedback loop. In this way, the activity of 32 is regulated through a feedback loop that involves the competition of 32 and unfolded proteins for binding with free DnaK/J chaperone pool. On the other hand, the regulation of 32 stability proceeds as follows. The value 32 is rapidly degraded (t1/2 = 1 min) during steady-state growth, but is stabilized for the first 5 min after temperature upshift. The chaperone DnaK and its co-chaperone DnaJ are required for the rapid degradation of 32 by the HS protease FtsH. RNAP bound 32 is protected from this degradation. Furthermore, FtsH itself being a product of the heat-shock protein expression, experiences a synthesis rate that is tied to the transcription/translation rate of DnaK/J. Therefore, as protein unfolding occurs, 32 is stabilized through a temporary relief of its sequestration from DnaK. However, as more proteins are refolded, and as the number of FtsH itself increases, there is a decrease in the concentration of 32 to a new steady-state concentration that is dictated by the balance between the temperature-dependent translation of the rpoH mRNA and the level of 32 activity modulated by the heat shock protein chaperones and proteases acting in a negative feedback fashion. In this way, the FtsH-mediated degradation of 32 is feedback-regulated. We refer to this as the FtsH degradation feedback loop. Some of the implications of this regulation on intrinsic noise attenuation will be addressed in this work. A biological block diagram illustrating the various regulation mechanisms at work in the heat shock response system is shown in Fig. 4.
A reduced deterministic model for the heat shock response
In previous work, we have developed a detailed model describing the dynamics of the heat shock response in the bacterium E. coli (17,25,26). Although the accurate quantitative modeling of the heat shock response involves a large set of nonlinear differential-algebraic equations, it can be demonstrated that the basic functional modules in this response along with their qualitative behavior can be captured by a simplified model whose components are briefly described in what follows. For ease of notation, we will denote the chaperones by D, the protease by F, the 32 by S, and the unfolded proteins by Uf. Furthermore, we refer to the total concentration of X by Xt and the free quantity of X by Xf. We will assume that D is produced at a rate kd and is degraded at rate d. The quantity D can reversibly bind to the -factor S or to the unfolded protein Uf. In turn, the level of free f influences the rate of transcription of protein D. The itself is produced at a temperature-dependent rate 
(T) and degraded through a regulated mechanism involving D and F. The proteins are folded through the interaction with D and unfolded with a temperature-dependent rate K(T). The chemical reactions describing this network are naturally divided into two categories—fast and slow. The fast reactions are assumed to be in equilibrium with respect to the slow reactions. In terms of differential equations (mass action), those reactions translate to
 | (2) |
where Pfolded is the concentration of folded proteins, U:D is the complex formed by the binding of Uf to D, and S:D:F is the complex formed by the binding of S to F mediated by the binding to D. The total concentration of protein D, unfolded proteins Uf, F, and -S are related to their other forms (free and in complex with other molecules) through the mass balance equations
 | (3) |
where Pt is the total number of proteins in the cell, considered here to be constant. We eliminate U:D and S:D:F from Eq. 2 by utilizing the fact that the binding reactions are fast compared to protein production and degradation, and write algebraic expressions:
 | (4) |
Straightforward algebraic manipulations of Eqs. 3 and 4 yield the relationship for Sf of
Similarly, Df can be expressed in terms of Dt. The exact expression is quadratic. However, an approximation is possible if one notices that usually St << Dt (St = 25–40 molecules/cell while Dt 
104 molecules/cell). In this case, one can drop S:D from mass-balance Eq. 3, and the expression for Df simplifies to
and consequently
Further simplifications can be carried by observing that the dynamics of Uf are much faster that those for St and Dt and therefore, can also be assumed to be in quasi-steady state. Under these assumptions, the expressions in Eq. 2 become
 | (5) |
where we have defined 
and 
The value Pt is the total number of proteins in the cell (assumed to be constant), Ks, Kf, Ku are binding rate constants, and 0, d, s are degradation constants. The parameter values are given in Table 2. This model reproduces qualitatively the dynamics of the heat shock response system for a plausible set of biochemical parameters while containing all the major feedback loops involved in heat shock regulation. Specifically, it possesses the FtsH degradation feedback loop, whose action appears in the last term of the differential equation for St Eq. 5. Plots showing the time histories of 32, chaperones, and unfolded proteins in the heat shock system are shown in Fig. 5. These plots show the characteristic behavior of a wild-type heat shock response. The concentration of 32 sharply peaks shortly after the onset of heat, then declines as a result of the feedback regulation to a new steady-state value at the high temperature. The level of chaperones increases accordingly, resulting in refolding of proteins and the recovery of their level to a low value. It is worth noting here that even though the description of the reduced-order model seems to be based on empirical interactions, the same description could be exactly derived from the full model through various algebraic manipulations and approximations. We omit the details of such a derivation here. The model in Eq. 5 does not aim at representing the full, current view of the molecular mechanisms involved in the heat shock response in bacteria, which is known to involve a large number of interacting proteins and various levels of redundancy. It is rather intended to provide a qualitative account of the core regulatory levels in the heat shock system, and as such, is particularly well suited for both stochastic simulations and analytical manipulations.
View larger version (8K):
[in this window]
[in a new window]
|
FIGURE 5 Deterministic trajectories of the simplified heat shock model. (a) Number of 32 per cell. (b) Number of Dt per cell. (c) Number of unfolded proteins per cell. Heat shock occurs at time 200.
|
|
A stochastic model for the heat shock response
A detailed stochastic representation of heat shock system would involve the identification of the molecular reactions appearing in every term of Eqs. 2 and 4. These elementary biochemical reactions, along with the propensity of their occurrence, are given in Table 3.
We first verify that the behavior predicted by the deterministic model is recovered when using the SSA to simulate the stochastic model of Table 3. Specifically, we investigate whether the mean of the stochastic model corresponds to the solution of the deterministic rate equations. This is crucial in our analysis since we use the deterministic model to define elementary complex reactions and subsequently identify the deterministic model trajectory with that of the stochastic mean. A sample realization obtained for the detailed elementary reactions model using Gillespie's SSA is shown in Fig. 6 for the number of 32 and total chaperone Dt at low temperature. On the same plot, we show the low temperature deterministic trajectory obtained from integration of Eq. 5. This trajectory converges to a stable steady state which can be shown to be the only biologically plausible (real, nonnegative, and corresponding to measured levels of Dt and St in the cell) fixed-point of Eq. 5 for the set of parameters given in Table 2. The system has also one stable fixed-point at high temperature. It can be identified in Fig. 6.
It can be seen in Fig. 6 that for the chosen values of the parameters, the deterministic and stochastic approaches virtually coincide and that the only difference between the deterministic and stochastic time courses is the presence of random fluctuations in the latter. We have used different initial conditions, and in all the cases we have observed the behavior of the long-term solution is the same.
A reduced stochastic model of the heat shock response
An alternative description of stochastic biochemical kinetics can be based on elementary-complex reactions, as opposed to the detailed elementary reaction description (27).
The use of elementary-complex reactions is often motivated by the need to circumvent the stiffness that results from different timescales in a system, which often makes the SSA prohibitively slow. In the absence of systematic model reduction methodologies in the stochastic setting, Michaelis-Menten and various other well-established approximations are carried out for the deterministic system, with the resulting expressions defined as new propensities for the elementary reactions governing the remaining states. Quasi-steady-state type of approximations applied to the probability functions of the fast species have also been attempted (28). The general range of validity and connections between these two approaches are at best unclear and still are open research questions. Therefore, the use of such approximations should be justified on a case-by-case basis. In the case of the heat shock system, it is warranted by the sharp separation in time- and concentration-scales between the state variables. Furthermore, as our aim is to gain qualitative insight into the stochastic properties of the system and characterize analytically its noise-rejection properties and the role of the feedback loops in this noise rejection, the approach based on the elementary-complex reactions is an appropriate first step. To monitor its validity and the accuracy of its predictions, we check against the detailed elementary reactions model simulated using the SSA.
In the context of the heat shock response, compact kinetic expressions obtained after application of various approximations (e.g., quasi-steady-state) and appearing in the kinetic expressions in Eq. 5 are the adopted elementary complex reactions. The effect of transitions brought about by these reactions is given in Table 4. A typical trajectory corresponding to the simulation of these complex reactions is shown in Fig. 6 (blue line). The figure indicates that the use of this description with our choice of parameters generates a reasonable approximation and serves as the starting point for analytical investigation.
|
REGULATED VERSUS CONSTITUTIVE DEGRADATION IN THE HEAT SHOCK RESPONSE SYSTEM
|
|---|
In the wild-type heat shock system, FtsH-dependent degradation of the factor, 32, implements a negative feedback loop. We demonstrate here that this regulated degradation possesses inherent benefits in terms of biochemical noise rejection. To assess these benefits, we compare the wild-type model of the heat shock to a virtual mutant model in which the dynamic degradation of 32 is replaced by a static (constitutive) degradation term.
We again compare these two cases both numerically and analytically.
Exact numerical analysis using the SSA
Using the SSA and the detailed elementary reactions given in Table 3, we investigate the exact statistical properties of the wild-type heat shock system and the alternative scheme of constitutive degradation described above. To achieve constitutive degradation, we decouple the production of FtsH from 32 (hence breaking the feedback loop) and adjust the elementary reactions in Table 3 accordingly. We focus on the chaperone (Dt) level, since it is a crucial indicator of protein folding. Wide stochastic excursions in this level are not desirable and cause performance degradation under both normal and heat shock conditions. A histogram of Dt compiled for a sample of 2400 SSA realizations of the wild-type heat shock at low temperature is shown in Fig. 7 a. The histogram indicates that the chaperones have a close-to-symmetric distribution around the mean, which in turn closely coincides with the deterministic steady state. Fig. 7 b also shows the histogram resulting from an equal number of samples obtained from the detailed stochastic simulation of the mutant system. These histograms clearly indicate that the distribution of chaperones is indeed narrower in the regulated degradation than in the unregulated degradation case. As a sensible quantitative measure for the effect of feedback, we computed the relative change in the coefficient of variation (standard deviation divided by the mean), which in this case is the same as the relative change in the standard deviation in the system. In mathematical terms, this improvement is expressed as
 | (6) |
where r is the coefficient of variation for the regulated degradation system and c is the coefficient of variation for the constitutive degradation system. For the data in Fig. 7 c, ß 
30.4%, which indicates a 30.4% improvement achieved by feedback relative to constitutive degradation.
Analysis using the LNA
In this section, we use the LNA to provide analytical insight into the factors contributing to the superior noise-rejection behavior of regulated degradation. We point here to the fact that we do not expect exact quantitative agreement between the LNA and SSA results. Error in the LNA estimates can result from two sources: the approximations involved in the elementary-complex description, in addition to the approximations inherent in the LNA method itself. However, the insight provided by the LNA will prove valuable in identifying the importance of the regulated stability of 32 in attenuating intrinsic noise in the heat shock system. The hallmark of this noise attenuation mechanism will be characterized by an enhanced feedback strength, which in turn translates to lower value for the chaperone variance. Again, as in the simple example, the use of LNA correctly identifies this effect qualitatively.
The wild-type heat shock system
We now use the LNA to compute the covariance matrix C of the chaperone level. We use the elementary-complex reactions in Table 4 and form the stoichiometry matrix S as
and the propensity vector as
 | (7) |
The stationary Jacobian matrix As (1) for the wild-type is
 | (8) |
where
while the diffusion matrix Ds is diagonal and given by
where
Solving the Lyapunov Eq. 1, we obtain the stationary covariance matrix. The steady-state variance for Dt can be shown to be
 | (9) |
The degradation mutant heat shock system
To assess the benefits of regulated degradation, we compare the wild-type model of the heat shock (expressions in Eq. 5) to a virtual mutant model in which the dynamic degradation of 32 is replaced by a static (constitutive) degradation term. This can be achieved if, for example, 32 is degraded by the action of a protease (still FtsH) whose production is not regulated by 32, i.e., Ft = 
= constant >> St. Under the same assumptions governing the wild-type model in Eq. 5, the mutant model equations will be
 | (10) |
As in the wild-type model, the elementary-complex reactions for the mutant model can be identified. They are listed in Table 5.
To compute the statistics of the chaperone level using the LNA, we repeat the same exercise as in the previous section. We break the mutant systems in Eq. 10 into complex elementary reactions, and form the stoichiometry matrix S as
and the propensity vector for the mutant system as
 | (11) |
To establish a sound base for comparison, we require that the steady-state values for the system variables be equal in both wild-type and mutant, i.e., 
and 
We can easily achieve this by setting the rate of constitutive synthesis of the protease in the mutant to be equal to that resulting from the dynamic synthesis, i.e., 
The stationary Jacobian matrix Acs of Eq. 1 for the mutant is
 | (12) |
where
Notice that with 
 | (13) |
Also, the stationary diffusion matrices Dcs in the mutant system is equal to its equivalent, Ds, in the wild-type system for equal stationary values of the states. Solving the corresponding Lyapunov equation for the mutant system, we get the following expression for the stationary variance of the static degradation mutant
At the same time, using the relations of Eq. 13 in Eq. 9 yields for the dynamic degradation wild-type
 | (14) |
It is now a straightforward algebraic manipulation to verify that 
independent of the choice of parameters, therefore demonstrating the clear benefit of the dynamic regulation of 32 stability in providing a more tightly distributed chaperone level.
|
DISCUSSION
|
|---|
Engineering systems use negative feedback as a preferred noise-attenuating mechanism. Not surprisingly, gene regulatory networks that must also function reliably in the presence of noise, ubiquitously use feedback. In terms of its filtering properties, a simple negative feedback loop can act as a lowpass filter. More sophisticated types of feedback, such as integral feedback, can shape bandpass filters. Integral feedback is, for example, implemented in bacterial chemotaxis (29). Mechanistically, negative feedback can be implemented through different sets of molecular interactions.
For example, protein products can negatively autoregulate their own synthesis, activity, or stability. Each of these strategies can have its own unique properties in terms of the speed and efficiency of the response it implements, as well as the stochastic properties of this response. The objective of this work was to offer a first glance at the stochastic properties of regulated degradation in the context of the heat shock response. We compared the wild-type model of the heat shock to a virtual mutant model in which the dynamic degradation of 32 is replaced by a static (constitutive) degradation term. We argued that this can be achieved if, for example, 32 is degraded by the action of a protease (still FtsH), which is not part of the 32 regulon. The mutant model where regulated degradation is impaired showed inferior performance in the presence of biochemical noise. This translated to chaperone levels that showed wide stochastic fluctuations around their desired mean value. On the other hand, the wild-type heat shock with intact 32 degradation showed smaller excursions around the same mean. We demonstrated that this behavior is the outcome of a feedback gain which assumed a higher value than in a constitutively degraded 32 mutant. This in turn resulted in lower variances of the protein levels in the presence of both intrinsic and extrinsic sources of noise. The power of this result is that it is demonstrated in a naturally occurring complex system such as the heat shock, and is therefore likely to hold true in similar systems or systems that share the same logic with the heat shock network.
Noise rejection aside, the regulated stability of proteins is also instrumental in implementing fast response transients. In a previous work, we have demonstrated this property in the heat shock system (17). We have shown that by promptly halting the degradation machinery, an immediate stabilization of 32 can be achieved. Coupled with an increase in the translation rate of 32, this strategy accounts for the sharp triggering and fast adaptation of the response.
Experimental and analytical verification of this prediction in addition to comparison of the two schemes in biologically relevant examples still await further investigation. Such investigation is very important since protein proteolysis is increasingly being characterized as a tightly regulated, specific, and temporally controlled process that plays a major role in a variety of basic cellular pathways. Examples are numerous and can be found across many organisms. In prokaryotes, the heat shock response is one of many examples where the control of protein degradation is an essential regulation strategy. Regulated proteolysis also abounds in eukaryotes. One notable example is the selective and programmed degradation of cell-cycle regulatory proteins such as cyclins, inhibitors of cyclin-dependent kinases, and anaphase inhibitors. This degradation control is essential for the timely progression of the cell cycle, which is driven by oscillations in the activity of cyclin-dependent kinases (Cdks). The activity of the Cdks is controlled by the intricate interplay resulting from the periodic synthesis and degradation of the positive regulators, cyclins, and the negative regulators, Cdk inhibitors. The different cyclins, specific for the G1, S, or M phases of the cell cycle accumulate and activate the Cdks at the appropriate times during the cell cycle. They are then degraded resulting in Cdk inactivation. Furthermore, levels of some Cdk inhibitors, which inhibit certain cyclin/Cdk complexes, also rise and decline at appropriate times again under appropriate degradation control by the ubiquitin system (9). Another example is the controlled ubiquitin-mediated degradation of tumor suppressors and proto-oncongenes, which forms the basis of healthy cell growth and proliferation. One of the best studied examples is in human cells, where the tumor suppressor p53 transcriptionally activates Mdm2, which in turn targets p53 for degradation (15). This mechanism is essential in reacting to stresses such as DNA damage. Indeed, after DNA damage, the activity of Mdm2 and its interaction with p53 decreases, therefore leading to stabilization of p53. The increase in p53 protein levels and in the transcription activity of p53 lead, in turn, to increased production of Mdm2, which dampens the response in a negative feedback fashion, much like the interaction between 32 and its proteases in the heat shock response system. The proper operation of the p53–Mdm2 degradation loop is essential as disruption in this degradation causes pathological conditions, including malignant transformations (30). A third example is the controlled proteolysis of many transcription factors such as the nuclear factor B (NF–B), a ubiquitous inducible transcription factor involved in the central immune system, inflammatory, stress and developmental processes (16). The most studied of these factors is NF–B1, a heterodimer which is retained in a latent form in the cytoplasm of nonstimulated cells through its association with inhibitory molecules called inhibitors of B, IB. After exposure of the cell to a variety of extracellular stimuli such as stress and viral and bacterial products, IBs are degraded rapidly, thereby activating the NF–B1 heterodimer, which is then translocated to the nucleus where it resumes its transcriptional activity.
These instances offer strong evidence that regulated protein proteolysis is an abundant and ubiquitously used strategy in gene regulatory networks. Consequently, experimental and mathematical investigation of its biological relevance should proceed hand-in-hand, to provide a full functional characterization of this important genetic regulation strategy.
|
ACKNOWLEDGEMENTS
|
|---|
This material is based upon work supported by the National Science Foundation under grant No. NSF-ITR CCF-0326576.
Submitted on February 1, 2005;
accepted for publication December 8, 2005.
|
REFERENCES
|
|---|
1. Thattai, M., and A. van Oudenaarden. 2001. Intrinsic noise in gene regulatory networks. Proc. Natl. Acad. Sci. USA. 98:8614–8619.[Abstract/Free Full Text]2. Paulsson, J. 2004. Summing up the noise in gene networks. Nature. 427:415–418.[CrossRef][Medline]
3. Morishita, Y., and K. Aihara. 2004. Noise reduction through interaction in gene expression and biochemical reaction processes. J. Theor. Biol. 228:315–325.[CrossRef][Medline]
4. Kepler, T., and T. Elston. 2001. Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations. Biophys. J. 81:3116–3136.[Abstract/Free Full Text]
5. Das, N., M. Valjavec-Gratian, A. Basuray, R. Fekete, P. Papp, J. Paulsson, and D. Chattoraj. 2005. Multiple homeostatic mechanisms in the control of p1 plasmid replication. Proc. Natl. Acad. Sci. USA. 102:2856–2861.[Abstract/Free Full Text]
6. Paulsson, J., K. Nordstrm, and M. Ehrenberg. 1998. Requirements for rapid plasmid cole1 copy number adjustments: a mathematical model of inhibition modes and RNA turnover rates. Plasmid. 39:215–234.[CrossRef][Medline]
7. Paulsson, J., and M. Ehrenberg. 2000. Random signal fluctuations can reduce random fluctuations in regulated components of chemical regulatory networks. Phys. Rev. Lett. 84:5447–5450.[CrossRef][Medline]
8. Becskei, A., and L. Serrano. 2000. Engineering stability in gene networks by autoregulation. Nature. 405:590–593.[CrossRef][Medline]
9. Coepp, D. M., J. Harper, and S. Elledge. 1999. How the cyclin became a cyclin: regulate proteolysis in the cell cycle. Cell. 97:431–434.[CrossRef][Medline]
10. Hilt, W., and D. Wolf. 1996. Proteasomes: destruction as a program. Trends Biochem. Sci. 21:96–102.[CrossRef][Medline]
11. Swain, P. 2004. Efficient attenuation of stochasticity in gene expression through post-transcriptional control. J. Mol. Biol. 344:965–976.[CrossRef][Medline]
12. Tsui, H., G. Feng, and M. Winkler. 1997. Negative regulation of mutS and mutH repair gene expression by the Hfq and RpoS global regulations of Escherichia coli K-12. J. Bacteriol. 179:7476–7487.[Abstract/Free Full Text]
13. Claverie-Martin, F., M. Diaz-Torres, S. Yancey, and S. Kushner. 1991. Analysis of the altered mRNA stability gene from Escherichia coli. J. Biol. Chem. 266:2843–2851.[Abstract/Free Full Text]
14. Varshavsky, A. 1997. The ubiquitin system. Trends Biochem. Sci. 22:383–387.[CrossRef][Medline]
15. Ryan, K., A. Phillips, and K. Vousden. 2001. Regulation and function of the o53 tumor suppressor protein. Curr. Opin. Cell Biol. 13:332–337.[CrossRef][Medline]
16. Baeuerle, P., and D. Baltimore. 1996. Nf-b: ten years after. Cell. 87:13–20.[CrossRef][Medline]
17. El-Samad, H., H. Kurata, J. Doyle, C. Gross, and M. Khammash. 2005. Surviving heat shock: control strategies for robustness and performance. Proc. Natl. Acad. Sci. USA. 102:2736–2741.[Abstract/Free Full Text]
18. van Kampen, N. 1992. Stochastic Processes in Physics and Chemistry. Elsevier Science Publishing, Dordrecht, The Netherlands.
19. Gillespie, D. T. 1992. Markov Processes: An Introduction for Physical Scientists. Academic Press, New York.
20. Gillespie, D. 1992. A rigorous derivation of the chemical master equation. Physica A. 188:404–425.[CrossRef]
21. Gillespie, D. 1976. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22:403–434.[CrossRef]
22. Elf, J., and M. Ehrenberg. 2003. Fast evaluation of fluctuations in biochemical networks with the linear noise approximation. Genome Res. 13:2475–2484.[Abstract/Free Full Text]
23. Tomioka, R., H. Kimura, T. Kobayashi, and K. Aihara. 2004. Multivariate analysis of noise in genetic regulatory networks. J. Theor. Biol. 229:501–521.[CrossRef][Medline]
24. Straus, D., W. Walter, and C. Gross. 1989. The activity of 32 is reduced under conditions of excess heat shock protein production in Escherichia coli. Genes Dev. 3:2003–2010.[Abstract/Free Full Text]
25. El-Samad, H., M. Khammash, H. Kurata, and J. Doyle. 2002. Robustness analysis of the heat shock response in E. coli. In Proceedings of the American Control Conference. Anchorage, AK. 1742–1747.
26. Kurata, H., H. El-Samad, T. Yi, M. Khammash, and J. Doyle. 2001. Feedback regulation of the heat shock response in E. coli. In Proceedings of the 40th IEEE Conference on Decision and Control. Orlando, Fl. 837–842.
27. Keizer, J. 1987. Statistical Thermodynamics of Nonequilibrium Proces
TOP
ABSTRACT
INTRODUCTION
) and regulated synthesis (
) models. Histograms were constructed using 820,000 samples. The parameter values in panel a are ß = 100, 
