INTRODUCTION

TOP ABSTRACT INTRODUCTION SINGLE GENE, NO FEEDBACK REGULATED SYSTEMS I: SELF-... BIFURCATIONS ESCAPE TIMES REGULATED SYSTEMS II: MUTUAL... DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D REFERENCES

The rate at which proteins are synthesized from individual genes is tightly regulated. In prokaryotes, this regulation is accomplished in part by the binding of regulatory proteins to stretches of DNA upstream (by definition) of the protein-coding region of the gene. Regulatory proteins can either inhibit or facilitate the binding of RNA polymerase to DNA or facilitate the isomerization of the DNA-RNA polymerase complex into a transcriptionally competent state. RNA polymerase processes along the DNA, transcribing the DNA into messenger RNA (mRNA). A ribosome can associate with mRNA and begin translation of the mRNA into an amino acid sequence as soon as a complete ribosome binding site emerges from the RNA polymerase. Translation can occur many times per transcript. The canonical introduction to this subject remains the monograph of Ptashne (1992).

In eukaryotes, regulatory proteins, referred to as transcription factors, are also used as one method of controlling gene expression. Another form of eukaryotic regulation occurs through chromatin-DNA interactions. In general, transcriptional regulation in eukaryotes is more complicated than in prokaryotes. However, we believe the results presented in this manuscript are relevant to both types of cells.

In this manuscript, we adopt the terminology used in studying prokaryotes and use "operator" to refer to upstream regulatory DNA sites. The term "promoter" refers to the nucleotide sequence to which RNA polymerase binds to begin transcription. Single genes may have multiple operators that can overlap with the promoter. The operator is said to be in an occupied state if a regulatory protein is bound to it and unoccupied otherwise. Chemical reactions that change the state of the operator are referred to as operator fluctuations. One of the main goals of this manuscript is to understand the role of operator fluctuations in transcriptional regulation. The mathematical methods developed here also directly apply to eukaryotic systems that regulate gene expression through transcription factors or chromatin-DNA interactions.

Much previous mathematical modeling of transcriptional regulation represents the production of gene product as a deterministic process (for a review, see Hasty et al., 2001b). In these models, both the gene-product concentration and operator state are treated as continuous. The probability of an occupied operator is given by a function of the regulatory protein concentration that is computed using thermodynamic arguments (Shea and Ackers, 1985).

In fact, however, there is now considerable experimental evidence that indicates the presence of significant stochasticity in transcriptional regulation in both eukaryotes and prokaryotes. In particular, several recent experiments on mammalian cells have supported the idea that gene initiation in response to an inductive signal is a stochastic process (Weintraub, 1988; van Roon et al., 1989; Fiering et al., 1990; Ko et al., 1990; Dingemanse et al., 1994; Walters et al., 1995). Additionally, there is evidence that chromatin-regulated gene expression is stochastic (Ahmad and Henikoff, 2001; Wijgerde et al., 1995) and that the initiation and deactivation of pigment expression during melanocyte differentiation proceeds in a random fashion (Bennett, 1983). Finally, studies on engineered gene circuits, which have been designed to act as toggle switches and oscillators, have revealed large stochastic effects (Elowitz and Leibler, 2000; Gardner et al., 2000; Becskei et al. 2001).

Several recent papers have reported theoretical investigations into the effects of fluctuations in gene regulation. Thattai and van Oudenaarden (2001) used simple models of transcription and translation to derive expressions for the means and variances of mRNA and protein number that compare favorably with experimental results. Using Monte Carlo simulations, McAdams and Arkin (1997) studied a more detailed model of these processes that also takes into account ribosome-RNase binding competition. Neither of these two investigations considered fluctuations in the state of the operator. Simple models that take into account fluctuations between active and inactive genetic states have been studied and used to explain induction-level heterogeneity among individual cells expressing steroid-inducible genes (Ko, 1991, 1992), and to support the conjecture that haploinsufficiency disease arises from stochastic gene expression (Cook et al., 1998). In these models, transcription and translation proceed deterministically while the gene is on. Additionally, these investigations relied on computer simulation and did not fully explore the consequences of fluctuations in regulated systems.

Numerical simulations of complete genetic networks have also been carried out. For instance, Arkin et al. (1998) considered a detailed stochastic model for the initial decision between the two developmental pathways (lysis and lysogeny) of bacteriophage . However, in this investigation the chemical kinetics of the operator fluctuations was assumed to be fast. This assumption allowed the operator states to be treated deterministically using a quasi-steady-state approximation. The role of noise has also been considered in engineered gene networks (Hasty et al., 2000, 2001b). In this work, fluctuations were added post-hoc to deterministic rate equations, and, therefore, the noise strength was an adjustable parameter.

Our aim here is to consider the origins of intrinsic fluctuations and to develop appropriate representations for them within the context of de novo stochastic models. We further consider a set of simplifying approximations found in various parameter limits, eventually arriving at deterministic models. Because our method starts from a microscopic description of the process, all the parameters in the approximate models are defined in terms of the underlying chemical rate constants. Our approximate schemes are important for two reasons. First, they provide insight into the dynamics of the system that cannot be gained from Monte Carlo simulations of the full model; and, second, numerical simulations using the approximate schemes can run orders of magnitude faster than Monte Carlo simulations of the full process.

In the first two sections and in the appendices, we present the mathematical techniques for manipulating the models. These methods are then used to examine the consequences of fluctuations in regulatory systems: 1) Noise destabilizes genetic switches. We compute the mean first passage times for a simple single-gene switch and examine their behavior under variations of several parameters, particularly the rate of operator transitions. 2) The destabilization of switches makes it necessary to generalize the notion of bifurcations. We examine stochastic bifurcations (Horsthemke and Lefever, 1984), qualitative changes in the probability density function under changes in parameters, again, with special attention to those bifurcations induced solely by changes in the rate of operator transitions.

We start by considering a simple model without feedback to establish our methods, then move on to a switching system consisting of a single self-promoting gene, and finally to a switching system composed of two mutually repressing genes.

	SINGLE GENE, NO FEEDBACK

TOP ABSTRACT INTRODUCTION SINGLE GENE, NO FEEDBACK REGULATED SYSTEMS I: SELF-... BIFURCATIONS ESCAPE TIMES REGULATED SYSTEMS II: MUTUAL... DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D REFERENCES

We begin with a model for a gene that has no feedback, direct or indirect, onto its own transcriptional regulation. Similar models have been used to study steroid-inducible genes (Ko, 1991, 1992) and haploinsufficiency disease (Cook et al., 1998). The simplicity of this model allows us to present the techniques we use for the analysis of regulated systems in a direct manner, uncomplicated by nonlinearities. The state space of the model consists of the number of gene product monomers (an integer variable) and the state of the gene's operator (binary). We use the term "gene product" to account for the lumping of mRNA and protein in this treatment.

The events that occur in this model can be represented as biochemical reactions (Fig. 1), and the master equation derived directly from that representation. In this manuscript, we adopt the following notational conventions. Uppercase calligraphic letters denote a molecule of a particular protein species or an operator (e.g.,  for a monomer and  for an operator). Uppercase letters represent state variables that denote the current number of molecules of a particular protein species or the current chemical state of the operator. Lowercase letters denote any allowed value of their uppercase counterparts. Because the state variables are pure numbers, all the rate constants have units of inverse time and are reported here as per second. Concentrations are formed by dividing the mean number of a protein species by the relevant volume (e.g., the cell volume). In this case, second-order rates constants would have units of 1/(time × concentration).

View larger version (18K): [in this window] [in a new window]   FIGURE 1   A schematic diagram illustrating transcriptional regulation without feedback. The operator has two possible states, occupied and unoccupied, and fluctuates between them. If the operator is empty, protein monomers  are produced at a rate 0; if the site is occupied, the production rate is 1. The operator does not affect the degradation of the protein product, which occurs with rate . We assume that the activator concentration is constant.

The degradation of gene product is written as

(1)

where  represents the monomer form of the expressed protein,  is the degradation rate, and  is used to denote a protein sink. We also use  to denote a protein source.

We write the spontaneous transitions between operator states, denoted ₀ (occupied) and ₁ (unoccupied), as

(2)

where the rate K sets the time scale for this reaction. The constants k₀ and k₁ are dimensionless and constrained to obey k₀ + k₁ = 1.

The reaction for the production of protein is written as

(3)

where rate _s (s = 0 or 1) depends on the chemical state of the operator. It is important to recognize that this simple reaction is an effective reaction representing a large number of component reactions together making up transcription, translation of mRNA into polypeptides and the folding of these polypeptides into proteins.

In this treatment, we distinguish two sources of stochasticity: operator fluctuations, and the combined action of transcription per se and degradation of protein product. In both cases, the variability results from the inherent discreteness of the state space and the randomness in the dwell time between discrete reactions. (Below, we will consider a further source of variability in the reaction in which gene product monomers interact to form dimers, although our discussion will serve primarily to indicate why these dimer fluctuations can be neglected.)

At any given time t, the state of the system described by reactions 1-3 is specified by the number of monomer proteins M(t) and the chemical state of the operator S(t). S(t) is equal to 0 if the operator is unoccupied and 1 if it is occupied. Thus, the various states of the system can be written as the ordered pair (m, s), where m is a non-negative integer and s is either 0 or 1. M(t) and S(t) are random variables, because the chemical reactions that change the state of the system occur randomly in time. If we assume that the dwell time in any particular chemical state of the system is exponentially distributed, then the system satisfies the Markov property. Physically, this means that the time evolution of the system is determined solely by its current state and is independent of its past. The Markov property allows us to write down a master equation for the time evolution of the probabilities p(t) = Pr[M(t) = m and S(t) = s] (van Kampen, 1992). Simply put, the master equation is a rate equation for p(t). Written out explicitly, the master equation for this process has the form

(4)

(5)

where we have suppressed the explicit time dependence of p(t) for readability. The first term on the right-hand side of Eqs. 4 and 5 represents the rate at which probability is flowing out of the state (m, s). It is the product of the average rate, (Kk_s + m + _s), at which transitions occur out of (m, s) and the probability p(t) of being in (m, s). Correspondingly, the other three terms on the right-hand side of these equations represent the rates at which probability flows into (m, s) from accessible states.

Although the notation used in Eqs. 4 and 5 makes these equations easy to interpret, it quickly becomes unmanageable as the systems become more complex. Therefore, we will adopt a more compact notation and combine these equations as

(6)

where the hatted index indicates "the other one":  = (s + 1) mod 1.

The partial moments, defined for integer j, are

(7)

Note that the zeroth moments m⁰_s are the marginal probabilities for the operator to be in state s. The moments are then the sums over operator states of the partial moments,

(8)

We can use Eq. 6 to derive ordinary differential equations (ODEs) for the time evolution of the partial moments (van Kampen, 1992). Because the reactions are all zeroth- and first-order, the ODEs for the partial moments factorize into independent pairs of linear equations, which can easily be solved. Here we are interested only in the steady-state values of the first moment and the variance,

(9)

and

(10)

where the overbar indicates the steady-state value. Under the given circumstances, we expect the first factor in the product on the right side of Eq. 10 to be of order one, the second factor to be of order m_*² and the order of the third to depend on the relative rates of product decay,  and operator transitions K. When product decay is much faster than operator transitions, the last factor is of order one and the variance is dominated by the second term, of order m_*². When the operator transitions are much faster, the last factor is very small, and the variance approaches the mean.

An appropriate measure of the relative size of the fluctuations is the coefficient of variation, which is defined as the ratio of the variance to the mean squared. The steady-state coefficient of variation can be written as

(11)

Note that, as the average number of monomers increases, the first term in the CV decreases. Even with large protein numbers, however, the CV can still be large due to fluctuations in the operator state. As these operator fluctuations become faster, the second term decreases as well. This stabilizing effect of fast operator fluctuations has been observed in computer simulations of stochastic gene expression (Cook et al., 1998). Below, we construct approximations to the dynamics for cases in which the product number is large or the operator fluctuations are fast.

Small-noise and fast-transition approximations to the master equation

For cases in which there is feedback regulation, we generalize Eq. 6 to include state-dependent transition rates and nonlinearities. When this is done, the moment equations no longer factorize. Numerical solutions can be difficult to obtain as well, especially when several genes are considered as part of a regulatory network. Therefore, we develop approximations to Eq. 6 that can be generalized immediately to the nonlinear cases of feedback regulation and are valid as one or the other source of variability becomes negligible.

Large steady-state gene-product level

(12)

(13)

(14)

u_x

(15)

(16)

(17)

(18)

Fast operator fluctuations

(19)

(20)

(21)

Simultaneous limits

(22)

(23)

(24)

(25)

(26)

(27)

(28)

	REGULATED SYSTEMS I: SELF-PROMOTER

TOP ABSTRACT INTRODUCTION SINGLE GENE, NO FEEDBACK REGULATED SYSTEMS I: SELF-... BIFURCATIONS ESCAPE TIMES REGULATED SYSTEMS II: MUTUAL... DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D REFERENCES

We now consider a system in which the gene product is itself an activating regulatory protein for its gene. This system is regulated by positive feedback. This type of regulation is thought to play a role in the developmental decision pathway of bacteriophage  (Ptashne, 1992) and has been used to construct a synthetic eukaryotic gene switch in Saccharomyces cerevisiae (Becskei et al., 2001). Here we study a minimal version of a system with positive feedback, which serves to illustrate the qualitative features of this type of regulation. For a more biologically complete treatment of the lysis/lysogeny decision pathway, the reader is referred to the pioneering work of Arkin et al. (1998). The system we consider is similar to that shown in Fig. 1, only now the activator is the gene's own product. After developing the models for this system and exploring the consequences of noise for bifurcations and spontaneous transitions in this system, we will come back to explore a second bipolar switch involving a pair of mutually repressing proteins.

Regulatory proteins often bind to their operators as dimers or higher-order oligomers (Ptashne, 1992). We assume that the active form of the protein in our system is a dimer. Thus, we explicitly consider the dimerization reaction and the stochasticity associated with it.

Letting  represent a protein dimer, the dimerization reaction is

(29)

where  is the dimensionless equilibrium dissociation constant, and  is the forward rate constant.

The operator transition reaction, Eq. 2, is now modified to explicitly include the role of the protein dimer

(30)

where K is now the forward transition rate and  is the dimensionless dissociation constant for this reaction. The other reactions remain as given above, except that we must specify the rate at which dimers are degraded. For simplicity, we have made the arbitrary choice that dimers are not degraded at all. That is, only the monomer form of the protein is unstable. Inclusion of an arbitrary dimer degradation is easily accommodated. Our choice does not have qualitative implications for this analysis within the range of parameter variation considered.

Let D(t) represent the number of dimers and N(t) the total protein number at time t. Then we have M(t) = N(t)  2D(t). The master equation for this process is

(31)

where p = Pr[N(t) = n, D(t) = d, and S(t) = s]. It is straightforward to generate sample paths of the "full" process described by the equations given above (see Fig. 6). These Monte Carlo simulations are usually computationally intensive, because the time scales of the various reactions involved can be very different. Therefore, we make use of the limiting cases discussed above to construct approximations to the master equation that are less computationally intensive and provide insight into the dynamics of the system.

The reaction given by Eq. 29 is generally assumed to be fast compared to all other reactions (McAdams and Arkin, 1997; Hasty et al., 2000, 2001a). Physically, this means that the monomer and dimer concentrations come to quasi-equilibrium before the total amount of protein changes appreciably. In Appendix A, we show how this assumption, along with the assumption that the coefficient of variation of D conditional on N is small, allows us to eliminate the dimer concentration from the problem. In this limit, we get the master equation, written in terms of the marginal density on M (In an abuse of notation, we let the marginal on M be represented as p and rely on the subscript itself to distinguish this probability distribution from that of N).

(32)

To arrive at this simple set of equations, it is necessary to assume that the dissociation constant for dimerization is large. Although this is not necessarily the regime of biological interest, it serves to simplify the presentation considerably and is the regime commonly used, usually implicitly, by other researchers. As discussed in detail in Appendix A, all the analyses presented below can be carried through when this assumption is relaxed; the qualitative nature of the results do not change. A direct comparison with experimental data, however, requires the more general treatment presented in Appendix A.

To move to the diffusion limit, we change to the dimensionless variables u = t and X = M/m_o, where m_o = ₁/ is the steady-state value of the mean monomer number for the system locked in the occupied state. In terms of these variables, the diffusion limit for large protein number produces

(33)

where the rescaled parameters are a₁ = 1, a₀ = ₀/₁,  = K/(³) and b = ²/.

As m_o   the fluctuations in the monomer concentration become negligible; Eq. 33 loses its diffusive terms and can be written as

(34)

All of the stochasticity in this system now derives from the operator fluctuations. The key difference between the above equation and Eq. 17 is the appearance of the x² factor multiplying ₀(x) in Eq. 34. As will be seen, this factor is responsible for the bistable behavior observed in the macroscopic limit of this system.

One useful feature of Eq. 34 is that an explicit expression for the steady-state marginal density  = ₀ + ₁ can be found,

(35)

where  is a normalization constant. Below, we discuss how this distribution undergoes bifurcations as a result of the fluctuations in the operator state.

Next, we consider the small-noise limit of the fluctuations in the monomer concentration and the fast-noise limit of fluctuations in the operator state. That is, both m_o and  are taken to be large, but finite. In Appendix B, we present a general algorithm for deriving an effective diffusion equation for the marginal density, and find

(36)

where

(37)

(38)

The steady-state solution to Eq. 36 is given by Eq. 25, using the above expressions for A and B.

Finally, in the deterministic limit , m_o  , we are left with the ODE

(39)

where we have introduced the potential

(40)

Loosely speaking, (x) can be thought of as an effective free energy for the system. Its local minima represent stable steady states of the concentration. The local maxima are energetic barriers that must be surmounted by thermal activation. In general, such a free energy function does not exist for nonequilibrium systems, as is the case for the mutual repressor system considered below.

	BIFURCATIONS

TOP ABSTRACT INTRODUCTION SINGLE GENE, NO FEEDBACK REGULATED SYSTEMS I: SELF-... BIFURCATIONS ESCAPE TIMES REGULATED SYSTEMS II: MUTUAL... DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D REFERENCES

The deterministic system given by Eq. 39 acts as a switch in the appropriate parameter regime: the system has two distinct stable fixed points (for a discussion of the necessary conditions required to make a biological switch, see Cherry and Adler, 2000). Figure 2 shows the bifurcation diagram for this system as a function of the two parameters b and a₀. The corresponding potentials (x) are also drawn on the diagram to illustrate the number of stable fixed points in each region. In the region where there are two stable fixed points, the system acts as a genetic switch.

View larger version (12K): [in this window] [in a new window]   FIGURE 2   The deterministic bifurcation diagram for the self-promoter. The dimensionless parameters are defined as b = 2/ and a0 = 0/1. The points 1, 3, and 4 indicate the parameter values used to produce the potentials shown in the figure. Points 1 and 2 specify the systems investigated in the Monte Carlo simulations shown in Fig. 4, both of which are monostable.

For stochastic systems, the notion of "stable fixed point" in the state space is not well defined. To generalize the idea of bifurcations for applicability to our stochastic models in such a way as to be consistent with the usual meaning as the fluctuations become vanishingly small, we focus our attention on the gain and loss of critical points in the probability density function (Horsthemke and Lefever, 1984). So, for example, a bifurcation in which a single stable fixed point becomes a pair of stable fixed points and a single unstable point corresponds to the transformation of a unimodal probability density function to one that is bimodal.

We now illustrate how fluctuations in the operator state change the dynamics of the system. In particular, these fluctuations can either induce bistability in regions that are deterministically (i.e., in the zero-noise limit) monostable or wash out regions of bistability. Figure 3 A is a bifurcation diagram for the stationary distribution given by Eq. 35 as a function of b and a₀. For this figure  = 50. The dashed curve corresponds to the deterministic bifurcation diagram shown in Fig. 2. The various regions of this graph are labeled with the numbers 1-6. In Fig. 3 B, qualitative features of the steady-state distributions for each region of Fig. 3 A are shown. As the vertical line shown in Fig. 3 A is crossed from left to right, an integrable singularity occurs at the lower boundary of the distribution. Likewise, as the horizontal line is crossed from top to bottom, an integrable singularity occurs at the upper boundary of the distribution. The points labeled 1 and 2 in Fig. 2 are monostable in the deterministic limit. In Fig. 3 B, however, we see that finite operator fluctuations induce a type of bistability. The stationary distribution has a local maximum at high concentrations and the distribution is singular at x = a₀. Note that there is small region near the cusp of the deterministic bifurcation curve where the deterministic system predicts bistability, but the fluctuations cause this behavior to be lost. Not surprisingly, the Monte Carlo simulations presented below reveal that fluctuations in the monomer concentration can also wash out bistable behavior. In region 3, the distribution is bimodal. This region grows and shifts as  is increased until it coincides with the deterministic limit (dashed curve). This effect can be seen in Fig. 3 C, which is a bifurcation diagram for the steady-state distribution as a function of b and  with a₀ = 0.05. This value of a₀ is used in the Monte Carlo simulations discussed next.

View larger version (21K): [in this window] [in a new window]   FIGURE 3   (A) The bifurcation diagram for Eq. 35 with  = 50. The regions numbered 1-6 correspond to qualitatively different steady-state distributions. The dashed curve is the bifurcation diagram in the deterministic limit shown in Fig. 2. (B) The steady-state distributions associated with the regions numbered in (A) and (C). (C) The bifurcation diagram for b versus  with a0 = 0.05.

Figure 4 shows the results of Monte Carlo simulations consistent with Eq. 32. The parameter values in the upper pair of panels correspond to the point labeled 1 in Fig. 2. In the deterministic limit, the system is monostable. From Fig. 3 A, however, we see that the finite operator fluctuations have induced bistability. This can be seen from the dashed curve in the top right panel of Fig. 4, representing the steady-state distribution when monomer fluctuations are ignored. In the left panel, a typical time series for the process is shown. There is no indication of bistability. The solid curve shown in the top right panel is the steady-state distribution when the small-and-fast noise approximations are used (i.e., Eqs. 36 and 25). By comparing the dashed and solid curves, we see that, even with this large value of m_o, monomer fluctuations can still be detected. Indeed, these fluctuations are responsible for washing out the bistable nature of the system. We find excellent agreement between this distribution and the histogram from the Monte Carlo simulations.

View larger version (39K): [in this window] [in a new window]   FIGURE 4   Sample paths and distributions for the self promoter. The dimer concentration has been eliminated using the quasiequilibrium approximation (Eq. 32). In all the panels,  = 10,000, 1 = 1000 s1, 0 = 50 s1, and  = 10 (a0 = 0.05). Each histogram corresponds to the time series on its left. Together, they illustrate the bifurcations that occur as  or K are varied. The bifurcation that occurs as K is varied is due solely to fluctuations in the operator state and does not occur in the macroscopic limit. The distributions shown as solid lines are the results of the small-and-fast noise approximation. The dashed curves are the steady-state distributions given by Eq. 35.

Comparison of the histograms reveals that bistability can be induced by varying either  or K. Experimentally,  could be altered, for example, by the addition of a protease. Comparing the middle panels to Fig. 3 C, we expect that this system will be bistable, and, indeed, the time series shown in the middle left panel reveals that this is the case, although low monomer levels are very unlikely. The lower state can be further stabilized by increasing . In the bottom panel, K has been reduced while  is unchanged from the top panel. The bistability of the system is evident. Note that, although both approximations do a pretty good job of reproducing the steady-state distributions, there is some discrepancy at low concentrations. This is exactly where we expect the approximations to break down.

Figure 5 is the same as Fig. 4 except that, instead of performing Monte Carlo simulations of the discrete process, we have generated sample paths based on the stochastic differential equation associated with Eq. 36. This method runs roughly an order of magnitude faster than the discrete Monte Carlo simulations. As can be seen, the agreement between the diffusion approximation and the full simulation is good, and it looks as if the diffusion approximation is faithfully capturing the dynamics as well as the steady-state distribution. We expand on this point in our discussion of the mean first passage time below.

View larger version (36K): [in this window] [in a new window]   FIGURE 5   Sample paths of the stochastic differential equation associated with the small-and-fast noise approximation and corresponding steady-state distributions. The solid curves are the same as in Fig. 4.

Finally, we illustrate the validity of ignoring the inherent dimer concentration fluctuations. Figure 6 A shows sample paths for M(t) and D(t) generated by Monte Carlo simulation of the process described by Eq. 31. In the deterministic limit, the system is bistable and this is indeed evident in the time series. The intrinsic fluctuations, however, allow for transitions between the high and low protein levels. Figure 6 B shows a comparison between the corresponding histogram and the steady-state distributions found using simultaneous application of the small-and-fast noise<