Stochastic modelling of nucleocytoplasmic oscillations of the transcription factor Msn2 in yeast

Didier Gonze, Michel Jacquet, Albert Goldbeter

Abstract

Stress induces oscillatory nucleocytoplasmic shuttling of the transcription factor Msn2 in yeast. The subcellular localization of Msn2 is controlled by the cAMP-dependent protein kinase, PKA. Recent experimental observations corroborated by a deterministic computational model for the cAMP–PKA pathway in yeast suggest that the oscillatory dynamics of Msn2 results from the periodic activation of PKA associated with stress-induced oscillations in the level of cAMP. The model accounts for the occurrence of oscillations of Msn2 in a window bounded by two critical values of the stress intensity. In contrast to the rather irregular oscillatory behaviour observed within single yeast cells by means of fluorescence measurements, the deterministic model can only produce a regular pattern of oscillations. To investigate whether the experimentally observed variability could be explained by molecular noise due to the small number of molecules involved in the oscillatory mechanism, we examine a stochastic version of the model for periodic nucleocytoplasmic shuttling of Msn2 coupled to oscillations in the cAMP–PKA pathway. The results of stochastic simulations compare well to the irregular oscillations observed experimentally in the nucleocytoplasmic shuttling of Msn2 in individual yeast cells. The stochastic model retains the property of oscillations within a range bounded by two critical values of stress intensity. We determine the dynamic behaviour as a function of this control parameter and show that the effect of noise markedly depends on the distance from the bifurcation points in the domain of oscillatory behaviour. Finally, we assess the role played by thresholds due to zero-order ultrasensitivity in phosphorylation–dephosphorylation cycles, both in the cAMP–PKA pathway and in the reactions controlling nucleocytoplasmic shuttling of Msn2. In regard to these thresholds, stochastic simulations show that large-amplitude variations of Msn2 associated with large-amplitude oscillations in cAMP can occur outside the domain of sustained oscillations predicted by the deterministic approach.

Keywords:

1. Introduction

The budding yeast protein Msn2 is an inducible transcription factor, which is activated by a large number of stress conditions (Martinez-Pastor et al. 1996; Gasch et al. 2000) such as heat shock (Boy-Marcotte et al. 1999), hyperosmotic shock, oxidative shock (Hasan et al. 2002), weak acid exposure (Gorner et al. 1998), ethanol addition and also by changes in available carbon sources (Boy-Marcotte et al. 1998; Gorner et al. 2002). This activation results in a rapid translocation of Msn2 from cytoplasm to nucleus (Gorner et al. 1998). The activation and shuttling of Msn2 are controlled by phosphorylation–dephosphorylation of the nuclear localization and nuclear export signals (NLS and NES). This phosphorylation is catalyzed by the cAMP-dependent protein kinase, PKA (Gorner et al. 1998, 2002). Recent experiments have shown that the control of the NLS and NES of Msn2 give rise to oscillatory shuttling of the transcription factor between the cytoplasm and nucleus (Garmendia-Torres et al. 2007). This phenomenon represents a new mode of dynamic behaviour based on NLS and NES control within the cell. The cellular population of Msn2 is able to shuttle in a coordinative manner between the cytoplasm and nucleus upon stress, which can take a variety of forms, including those listed above, exposure to light emitted from the microscope (Jacquet et al. 2003), and carbon source limitation (Medvedik et al. 2007). Oscillations in the subcellular localization of Msn2 have been attributed to oscillations of cAMP within the cell (Garmendia-Torres et al. 2007). The oscillations can be initiated by stress-induced changes in the cAMP–PKA pathway and likely originate from the negative feedback exerted by PKA upon the accumulation of cAMP (Nikawa et al. 1987).

In a previous report, a deterministic model was proposed to account for the oscillatory shuttling of Msn2 (Garmendia-Torres et al. 2007). This model shows that, on the basis of known regulation of the cAMP–PKA pathway, oscillatory dynamics are likely to occur over a wide range of parameter conditions. In particular, the model shows that oscillations of cAMP and Msn2 occur in a window bounded by two critical values of the stress intensity. Obviously, the deterministic model can only produce a regular pattern of oscillations, in contrast to the observed oscillatory behaviour within the cells, which is rather irregular. Fluorescence measurements reflect what happens in individual cells. Therefore, the source of variability is intracellular instead of originating at the level of the cell population. The question thus arises as to whether the experimentally observed variability could be explained by molecular noise due to the small number of molecules involved in the oscillatory mechanism that underlies the periodic nucleocytoplasmic shuttling of Msn2. To address the role of molecular noise, the model governed by deterministic laws is no more appropriate, and we need to resort to a stochastic description (Gillespie 1977, 2007; Barkai & Leibler 2000; Gonze et al. 2002a).

Here we investigate the dynamic behaviour of a stochastic version of the deterministic model recently proposed for Msn2 oscillatory shuttling based on oscillations in the cAMP–PKA pathway. We compare the results of stochastic simulations with the irregular oscillations observed experimentally in the nucleocytoplasmic shuttling of Msn2. We first recall in §2 the main results obtained for the deterministic model for stress-induced oscillations in the cAMP–PKA pathway and in the associated nucleocytoplasmic shuttling of Msn2 in yeast (Garmendia-Torres et al. 2007). In §3 we present the stochastic version of the model. The effect of molecular noise on Msn2 oscillatory shuttling is considered in §4. In §§5 and 6 we examine the effect of stress intensity and zero-order ultrasensitivity on stochastic oscillations. Section 7 is devoted to the discussion of these results.

2. Deterministic model for Msn2 nucleocytoplasmic oscillations based on oscillations in the cAMP–PKA pathway

2.1 Presentation of the deterministic model

The cAMP–PKA pathway in yeast is regulated by a strong negative feedback loop (Nikawa et al. 1987). To study the possibility that this regulation gives rise to sustained oscillations of cAMP, we recently considered a computational model for the cAMP–PKA system and investigated whether it could produce oscillations over a wide range of parameter values that match at least some experimentally determined values. The model is shown in figure 1 and presented in more detail in the electronic supplementary material that accompanies our original publication (Garmendia-Torres et al. 2007).

Figure 1

Scheme of the model for Msn2 nucleocytoplasmic shuttling coupled with the cAMP–PKA pathway in yeast. The cAMP–PKA system involves GEF, RAS, GAP, adenylate cyclase (CYCL), phosphodiesterase (PDE), cAMP and PKA. R2C2 represents the inactive form of PKA, while C represents its free, active catalytic subunit. For clarity of the scheme, the forms R2C2 and C are shown only in the cytoplasm but are also present in the nucleus. Four forms of Msn2 are considered: cytosolic non-phosphorylated Msn2 (MC), cytosolic phosphorylated Msn2 (MCP), nuclear non-phosphorylated Msn2 (MN) and nuclear phosphorylated Msn2 (MNP). Stress is assumed to decrease the activity of GEF and to enhance the dephosphorylation of Msn2. Two negative feedback loops in the cAMP–PKA pathway involve the activation of GAP and of PDE through phosphorylation by the catalytic subunit C of PKA. Oscillations in cAMP and PKA resulting from this regulation underlie the oscillatory nucleocytoplasmic shuttling of Msn2 (see text and Garmendia-Torres et al. (2007) for further details).

In yeast, cAMP is produced by an adenylate cyclase activated by Ras proteins (Toda et al. 1985), which are themselves activated by association with GTP, catalyzed by the guanine exchange factor, GEF (Broek et al. 1987). Active Ras is reset to its GDP form by a GTPase-activating protein, GAP (Tanaka et al. 1990). Binding of cAMP to the inhibitory subunit of PKA results in the dissociation and concomitant activation of the two catalytic subunits (C) of PKA (Zaremberg & Moreno 1996). cAMP is hydrolyzed by the two phosphodiesterases Pde1 and Pde2 of low and high affinity, respectively (Uno et al. 1983). A strong negative feedback loop is exerted by PKA on cAMP synthesis, as evidenced by the thousand-fold increase in cAMP in cells retaining only a weak activity of PKA (Nikawa et al. 1987). In the model (see Garmendia-Torres et al. 2007), for the sake of simplicity, we consider as single entities the redundant functional units RAS, GEF, GAP, PDE and C. Negative feedback on cAMP accumulation takes the form of PKA activation of both the phosphodiesterase, PDE, and the GTPase-activating protein, GAP, to hold with the original finding that the fully deregulated level of cAMP found in mutants with weak PKA activity was observed only in cells containing constitutively active RAS, and in which the two phosphodiesterases were deleted (Nikawa et al. 1987).

Coupling of Msn2 to cAMP variations is mediated through phosphorylation by PKA, which promotes cytoplasmic localization of the transcription factor. Msn2 can be phosphorylated by PKA both in the cytoplasm and in the nucleus since in glucose-grown cells kinase is present in both compartments (Griffioen et al. 2000) and cAMP is small enough to diffuse rapidly to the nucleus. Although the NLS is controlled by four phosphorylation sites (Gorner et al. 1998), we considered only one phosphorylation–dephosphorylation event in our simplified model because this is sufficient to couple the shuttling of Msn2 to the oscillatory dynamics of the cAMP–PKA system. The effect of stress is likely manifold: regardless of its precise nature—be it chemical, light emitted from the microscope, or carbon source limitation—we assumed that it controls the cAMP–PKA pathway through the inactivation of GEF (Wang et al. 2004), and enhances the activity of phosphatases acting on Msn2, as can be deduced from the experimental results showing that Msn2 still responds to stress by entering the nucleus in the absence of any variation in cAMP (Garmendia-Torres et al. 2007).

2.2 Equations of the deterministic model

The time evolution of the deterministic model is governed by the following set of ordinary differential equations (for further details on the equations and on the definition of the parameters, see Garmendia-Torres et al. 2007):Embedded Image(2.1a)Embedded Image(2.1b)Embedded Image(2.1c)Embedded Image(2.1d)Embedded Image(2.1e)Embedded Image(2.1f)Embedded Image(2.1g)withEmbedded Image(2.2a)

Embedded Image(2.2b)

Embedded Image(2.2c)

In the above equations, GEFt, GAPt, RASt, CYCLt, PDEt and PKAt denote, respectively, the total concentrations of the proteins GEF, GAP, RAS, adenylate cyclase (CYCL), PDE, and PKA; cAMP denotes the cellular concentration of cAMP. The negative feedback regulation by PKA is introduced in equations (2.1b) and (2.1e) by assuming that the catalytic subunit C of PKA catalyzes the phosphorylation of the inactive GAP protein, GAPi, into the active form GAPa, and activates similarly the phosphodiesterase, PDE. We assume that, regardless of its precise nature, stress acts by reducing the maximum rate of the enzyme that inactivates GEF; the dimensionless parameter Str measures stress intensity.

To couple Msn2 shuttling to the variations of cAMP, we incorporate into the model the phosphorylation of Msn2 by PKA. Then the time evolution of the fractions of non-phosphorylated and phosphorylated Msn2 in the cytosol (MC and MCP) and in the nucleus (MN and MNP) is governed by the following kinetic equations:Embedded Image(2.3a)Embedded Image(2.3b)Embedded Image(2.3c)Embedded Image(2.3d)where kt1, kt2, kt3 and kt4 denote first-order apparent rate constants for transport of the various forms of Msn2 into and out of the nucleus. V9 and V11 denote, respectively, the maximum rates of phosphorylation of the cytosolic and nuclear forms of Msn2 by PKA. We assume for simplicity that PKA and cAMP are present in equal amounts in both the cytosolic and nuclear compartments, and that the activation of PKA by cAMP described by equations (2.1g), (2.2a)–(2.2c) also takes place in the nucleus. If unequal amounts of PKA in the two compartments were to be considered, then C should be multiplied by a parameter in equations (2.3b) and (2.3c), to reflect such an asymmetry. Similarly, different equations could be considered for the time evolution of cytosolic and nuclear cAMP, instead of describing them by a single equation, as done with equation (2.1f). The assumptions retained also simplify the stochastic treatment (see §3 below). The maximum rates of dephosphorylation of the cytosolic and nuclear forms of Msn2 are assumed to be proportional to stress, so that V10 and V12 are multiplied by parameter Str.

The total concentrations of Msn2 in the cytosol (Mcyto) and in the nucleus (Mnucl) are given byEmbedded Image(2.4a)

Embedded Image(2.4b)

The conservation equation MSN2t=Mcyto+Mnucl allows us to compute the concentration of one of the four forms of Msn2 from the three other concentrations and to relinquish the corresponding kinetic equation in equations (2.3a)–(2.3d).

Equations (2.1a)–(2.4b) above are the same as in our previous publication (Garmendia-Torres et al. 2007), the only difference being that here we converted molar fractions into concentrations for all the variables considered, by multiplying molar fractions by the total concentrations of the corresponding proteins. This transformation facilitates the link with the stochastic version of the model, since it is straightforward to go from concentrations to numbers of molecules.

2.3 Dynamic behaviour of the deterministic model as a function of stress intensity

Depending on stress intensity measured by parameter Str, the model predicts that an oscillatory regime separates different types of steady-state situations corresponding to low or high levels of cAMP and PKA activity, associated with distinct patterns of Msn2 subcellular localization (see Garmendia-Torres et al. (2007) and figure 2). At low stress values, the system reaches a stable steady state corresponding to a high level of cAMP (figure 2a), which leads to the predominantly cytoplasmic localization of Msn2 (figure 2b). Beyond a critical value of stress intensity, oscillations of cAMP occur spontaneously (figure 2a). These oscillations are accompanied by oscillations of all variables of the cAMP–PKA pathway, including PKA, and by the periodic shuttling of Msn2 between the cytoplasm and nucleus (see Garmendia-Torres et al. (2007) and figure 3a). Finally, above a second, higher critical value of stress intensity, cAMP reaches a stable, low steady-state level (figure 2a), while Msn2 ceases to oscillate and becomes predominantly located in the nucleus (figure 2b).

Figure 2

Bifurcation diagram showing the envelope of oscillations for (a) cAMP and (b) nuclear Msn2 as a function of stress intensity (Str) in the deterministic model (Garmendia-Torres et al. 2007). Parameter values are given in table 1. The values of Str marked A–E correspond to the five situations considered in figure 7. Dots representing data obtained by numerical integration of equations (2.1a)–(2.1g) and (2.3a)–(2.3d) correspond to stable periodic or steady-state behaviour. In (a), dashed lines obtained with the program Auto (Doedel 1981) correspond to unstable periodic or steady-state behaviour.

Figure 3

Comparison of stochastic and deterministic oscillations. The curves show the oscillations of RGTP, active adenylate cyclase (CYCL), cAMP, active PKA and nuclear Msn2 obtained for Str=1 in the (a) deterministic and (b) stochastic versions of the model. Stochastic simulations were performed for Ω=200. Other parameter values are the same as those given in the legend to table 1.

The model predicts that when entering the domain of oscillations, the shuttling of Msn2 between the nucleus and cytoplasm is at first such that Msn2 remains predominantly in the cytoplasm. Then, upon increasing stress intensity, Msn2 oscillates between the cytoplasm and nucleus. Finally, at still higher stress intensities, before leaving the oscillatory domain, Msn2 oscillates with reduced amplitude and remains predominantly in the nucleus (see Garmendia-Torres et al. (2007) and figure 2b).

3. Stochastic model for oscillatory nucleocytoplasmic shuttling of Msn2

To assess the effect of molecular noise on nucleocytoplasmic oscillations of Msn2, we develop a stochastic version of the deterministic model outlined in §2. To this end, we must decompose the scheme of figure 1 into a sequence of biochemical steps. Two possibilities exist: either we consider a fully detailed decomposition into elementary reaction steps, or we adopt a semi-detailed description. In a previous study of stochastic models for circadian rhythms, we referred to these two approaches as the developed and non-developed versions of a stochastic model. For example, when considering a Michaelis–Menten reaction, we can either decompose it into the reversible formation of the enzyme–substrate complex followed by its catalytic decomposition to form the product—this is the fully detailed treatment—or we may skip the detailed decomposition and treat the probability of synthesis of the reaction product as simply being proportional to the Michaelis–Menten rate function. We compared the two approaches and showed that they give similar results in the case of circadian rhythms (Gonze et al. 2002b, 2003). The advantage of the non-developed stochastic approach is to markedly reduce the number of variables as well as the number of adjustable parameters and, hence, to reduce the complexity of stochastic simulations. Here, for the sake of simplicity, we will focus on a non-developed stochastic version of the model for Msn2 oscillatory shuttling linked to the oscillations in the cAMP–PKA pathway in yeast.

The reaction steps of the non-developed stochastic model are listed in table 1. Steps 1–10 pertain to the activation and subsequent inactivation of GEF, GAP, Ras protein, adenylate cyclase and phosphodiesterase. Steps 11 and 12 refer to cAMP synthesis by adenylate cyclase and degradation by phosphodiesterase. Steps 13 and 14 pertain to reversible PKA activation by cAMP, in the cytosol and in the nucleus. Finally, steps 15–22 refer to phosphorylation–dephosphorylation and reversible translocation of Msn2 between the cytoplasm and nucleus. For each step the reaction propensity is indicated in the right column. The probability of occurrence pi of reaction step i is then calculated from the reaction propensity wi according to the expression Embedded Image.

View this table:
Table 1

Reaction steps considered in the stochastic model, with their associated probability of occurrence. (Parameter values (taken from and defined in Garmendia-Torres et al. 2007) are: Str=1, GEFt=4 μM, V1=1 μM min−1, V2=1 μM min−1, K1=K2=GEFt×0.05 μM, GAPt=1.5 μM, V3=3.5 μM min−1, V4=1.3 μM min−1, K3=K4=GAPt×0.01 μM, RASt=250 μM, V5=240 μM min−1, V6=600 μM min−1, K5=K6=RASt×0.001 μM, CYCLt=0.7 μM, ka=RASt×0.01 μM−1 min−1, ki=1 min−1, ks=CYCLt×4 min−1, PDEt=0.5 μM, V7=3.333 μM min−1, V8=1.5 μM min−1, K7=K8=PDEt×0.01 μM, kd=100 min−1, Kmd=20 μM, PKAt=0.3 μM, a=1 μM−2 min−1, r=1 μM−2 min−1, MSNt=1 μM, V9=3.333 min−1, V10=0.6 μM min−1, K9=K10=MSN2t×0.05 μM, V11=3.333  min−1, V12=2 μM min−1, K11=K12=MSN2t×0.05 μM, kt1=10 min−1, kt2=0.001 min−1, kt3=0.001 min−1, kt4=10 min−1. Parameter Ω, which is used to convert the concentration into the number of molecules, has the unit of a volume (10−18 l). Parameters Ki (i=1, …, 12) are Michaelis constants given in dimensionless form in Garmendia-Torres et al. (2007) and expressed here in concentration units through multiplication by the total concentration of the corresponding substrate. Likewise, owing to the change from molar fraction to concentration, parameters ka and ks given in Garmendia-Torres et al. (2007) in min−1 have been multiplied here, respectively, by the total concentrations of Ras (RASt) and adenylate cyclase (CYCLt).)

The stochastic version of the model represented by the 22 reaction steps listed in table 1 has been studied numerically by means of the exact stochastic algorithm proposed by Gillespie (1977, 2007). Varying the number of molecules in the numerical simulations can be achieved by changing parameter Ω in the expressions for the reaction propensities of the various reactions in table 1. Thus, decreasing the value of Ω will decrease the number of molecules involved in the cAMP–PKA pathway and in the translocation of Msn2.

4. Effect of molecular noise on Msn2 oscillations

We present in figure 3b the results of numerical simulations of the stochastic model of table 1, for the value Str=1. This value corresponds to point C in the middle of the oscillatory range in the bifurcation diagram established in figure 2a as a function of stress intensity for the deterministic version of the model. The oscillations are compared with the predictions of the deterministic model (figure 3a). Shown are the time evolution of RGTP (active RAS), the active form of adenylate cyclase, cAMP, the catalytic subunit of PKA, and the nuclear form of Msn2. We see that all variables oscillate, either regularly in the deterministic model or in a noisy fashion in the stochastic approach. The abrupt nature of oscillations in RGTP is due to the fact that for the parameter values considered the RAS module operates under conditions of zero-order ultrasensitivity (Goldbeter & Koshland 1981). In these conditions (see §6), a sharp threshold indeed exists for the activation of Ras, which acquires an all-or-none nature (Garmendia-Torres et al. 2007).

Shown in figure 4a are typical oscillations of nuclear (and, in antiphase, cytoplasmic) Msn2 observed in the experiments. Figure 4b represents the time course of nuclear and cytoplasmic Msn2 predicted by the deterministic model in the domain of sustained oscillations. The most conspicuous difference between the curves in figure 4a and 4b is the irregular nature of the oscillations observed in the experiments but not, of course, in the deterministic simulations. The stochastic simulations illustrated in figure 4c show that the stochastic approach allows us to recover this important aspect of physiological oscillatory behaviour. The resemblance between the oscillations observed experimentally and those obtained by stochastic simulations is striking. This indicates that the fluctuations observed in individual yeast cells may be attributed to molecular noise.

Figure 4

Comparison of nucleocytoplasmic oscillations of Msn2 observed in (a) the experiments, (b) the deterministic model and (c) the stochastic model. The experimental time series is reproduced from fig. 1C in Jacquet et al. (2003). The deterministic oscillations have been obtained by numerical integration of the differential equations (2.1a)–(2.4b) (see also Garmendia-Torres et al. 2007). The stochastic time series has been obtained by means of the algorithm of Gillespie (1977, 2007) for Ω=75, which value yields the numbers of molecules of the order of those observed in yeast (Garmendia-Torres et al. 2007). Other parameter values are the same as given in the legend to table 1. Shown in (a–c) are the fractions of nuclear Msn2 (red curve) or cytosolic Msn2 (blue curve).

To further characterize the effect of molecular noise on the oscillations in the cAMP–PKA pathway and in the associated nucleocytoplasmic shuttling of Msn2, we determined the effect of parameter Ω, which allows us to vary the number of molecules involved in these cellular processes. Shown in figure 5 are the oscillations in cAMP (left column) and nuclear Msn2 (centre column) for (a) the deterministic model and (bd) the stochastic model with Ω decreasing from 250 to 25. When Ω is large, e.g. 250, the stochastic oscillations are more regular and resemble those observed in the deterministic case. As Ω progressively decreases, the number of molecules diminishes and the effect of noise becomes more significant, until noise eventually obliterates periodic behaviour. In figure 5ad (right column) we plot the histogram of periods obtained as a function of Ω. In the deterministic case the period is by definition constant and the histogram reduces to a vertical line. As expected the histogram becomes wider as the value of Ω decreases. The total number of molecules of Msn2 present in the system is obtained by multiplying by Ω the constant, total concentration of Msn2, MSNt, which is taken as equal to 1 μM (see legend to table 1). The comparison of the curves in figure 5c,d indicates that noise progressively overcomes coherent rhythmicity as the total number of Msn2 molecules decreases from 75 to 25 upon going from Ω=75 to 25. Note, however, that the number of cAMP and PKA molecules also decreases when changing the value of Ω, which also contributes to enhancing the effect of noise.

Figure 5

Effect of molecular noise on the oscillations of cAMP and Msn2. Oscillations in cAMP and in nuclear Msn2 are shown in the left and centre columns. The right column shows the histogram of periods determined for oscillations of nuclear Msn2 over 4×103 min. In the stochastic simulations, the number of molecules is varied by changing parameter Ω. (a) Deterministic oscillations are compared with stochastic oscillations obtained for (b) Ω=250, (c) 75 and (d) 25. Noise progressively overcomes the oscillations as Ω decreases to a value of the order of 25. Other parameter values are the same as given in the legend to table 1. The periods, or, rather, the time intervals between successive peaks, were determined by computing the time difference between two successive passages of nuclear Msn2 through the value 0.5 upwards. To filter out fluctuations around that threshold, periods smaller than (arbitrarily) 2 min were discarded.

We may decrease specifically the number of Msn2 molecules with respect to the number of molecules of all other species (including cAMP and PKA) by selecting in table 1 a smaller value for Ω, for steps 15–22 only. Then we obtain more robust oscillations of Msn2 (figure 6) than the oscillations obtained for the same number of Msn2 molecules in figure 5d (also established for Ω=25). The oscillations in cAMP (figure 6a) are indeed as robust as those shown in figure 5b (also established for Ω=250) and less noisy than those obtained for Ω=25 in figure 5. This result is further corroborated by the more rapid decrease in the autocorrelation of the time evolution in the case of figure 5 compared with the case of figure 6 (data not shown). The fact that Msn2 shuttling is driven by robust oscillations in the cAMP–PKA pathway thus enhances the robustness of Msn2 oscillatory shuttling with respect to molecular noise.

Figure 6

Stochastic oscillations obtained for a large number of molecules in the cAMP/PKA system relative to the number of Msn2 molecules. Shown are the oscillations in (a) cAMP and (b) nuclear Msn2, and (c) the histogram of periods determined for oscillations of nuclear Msn2 over 4×103 min. These results have been obtained by setting in table 1 Ω=Ω1=250 in reaction steps 1–14, and Ω=Ω2=25 in steps 15–22. To keep the same qualitative oscillatory dynamics, parameters V9 and V11 in table 1 were multiplied by the ratio Ω2/Ω1=0.1 and were thus taken equal to 0.333 min−1. The results show that the robustness of nucleocytoplasmic Msn2 oscillations can be increased if the cAMP/PKA signal is more robust, which happens when the number of molecules involved in the cAMP/PKA system increases. Other parameter values are the same as given in the legend to table 1.

5. Effect of stress intensity in the stochastic model for Msn2 oscillatory shuttling

To determine the effect of stress intensity (Str) on oscillations in the stochastic model, it is useful first to return to the dynamic behaviour of the deterministic model as a function of this key control parameter. We showed in figure 2 the envelopes of oscillations of (a) cAMP and (b) nuclear Msn2 as a function of parameter Str. As recalled in §2.3, the curves indicate that sustained oscillations occur in a domain bounded by two critical values of stress intensity corresponding to bifurcation points. The vertical lines refer to five values of Str, which correspond to situations where the system is just below (Str=0.6) or just above (Str=0.75) the lower critical value, well in the middle of the oscillatory domain (Str=1 or 1.3), and to the right of the oscillatory domain (Str=2). Holding the value of Ω constant we performed stochastic simulations for these five values of Str to determine the effect of the distance from a bifurcation point.

The results of these stochastic simulations are shown in figure 7 where rows (a)–(e) correspond to the values of Str marked by the corresponding capital letter in figure 2a. The time evolution of cAMP and nuclear Msn2 is plotted in the left and centre columns, respectively, while the corresponding deterministic phase plane trajectory is shown in the right column, where the asymptotic steady or periodic state is represented by a white dot or white closed trajectory (corresponding to a stable limit cycle), respectively. Below the lower bifurcation point the dynamic behaviour corresponds to fluctuations around a stable steady state (figure 7a). Just above the lower critical value of Str small-amplitude oscillations occur around an unstable steady state and the noisy trajectory surrounds the deterministic limit cycle of small amplitude (figure 7b, established for Str=0.75). When the system is located well inside the oscillatory domain (Str=1), periodic behaviour is clearly discernible (figure 7c) and the phase plane trajectory takes the form of a high-amplitude, noisy limit cycle surrounding the deterministic limit cycle trajectory (white curve).

Figure 7

The effect of stress intensity on the stochastic oscillations of cAMP and Msn2. The curves show the oscillations in cAMP (left column) and in nuclear Msn2 (centre column) obtained in the stochastic model for increasing values of stress intensity indicated by the corresponding vertical lines in the diagrams in figure 2: (a) Str=0.6, (b) 0.75, (c) 1, (d) 1.3 and (e) 2. Other parameter values are the same as given in the legend to table 1. The phase plane trajectory associated with the noisy oscillations is shown in the right column. The white points or curves indicate (a,e) the stable steady state or (b–d) the limit cycle trajectory predicted by the deterministic model. The deterministic results were obtained by integration of equations (2.1a)–(2.3d) for the parameter values listed in table 1.

For a higher value of stress intensity (Str=1.3) the deterministic bifurcation diagram of figure 2a shows that the oscillations of cAMP have nearly the same amplitude as for Str=1, but the amplitude of Msn2 oscillations is somewhat reduced. This property is reflected by stochastic simulations, which indicate that the phase plane trajectory takes the form of a noisy limit cycle of smaller amplitude, restricted to the portion of the plane corresponding to a predominantly nuclear location of Msn2 (compare figure 7c and figure 7d). Finally, when the system operates to the right of the oscillatory domain, above the upper critical value of Str, the system settles again in a stable steady state and the effect of noise is similar to that observed to the left of the oscillatory domain, but the situation is inverted since Msn2 is now predominantly nuclear (figure 7e).

6. Effect of zero-order ultrasensitivity

When comparing the amplitude of stochastic changes of cAMP and nuclear Msn2 in figure 7, we see that except in (c) and (d) established in the middle of the oscillatory domain, large peaks of cAMP are often not reflected by large-amplitude changes in Msn2 subcellular localization. This is due to the fact that the amplitude of variations in PKA activity predicted by the deterministic and stochastic versions of the model is smaller than the amplitude of cAMP variations (see figure 3 in the case of oscillations for Str=1). The relatively reduced periodic changes in PKA are nevertheless capable of producing large-amplitude oscillations in nuclear Msn2 in figure 7c,d because a threshold exists in the control of Msn2 localization by PKA, owing to zero-order ultrasensitivity in Msn2 phosphorylation–dephosphorylation (Goldbeter & Koshland 1981; Garmendia-Torres et al. 2007).

If this explanation is correct, we expect that the amplitude of Msn2 oscillations should progressively decrease when zero-order ultrasensitivity becomes less and less present. This is indeed what we observe in figure 8 when the values of the Michaelis constants K9K12 increase (in μM) from (a) 0.1, to (b) 0.5 and (c) 5. The situation corresponding to the value 0.05 is shown in figure 3. Only when the Michaelis constants are smaller by at least one order of magnitude than the total concentration of Msn2 (e.g. in figure 3 or in figure 8a)—which conditions favour ultrasensitivity (Goldbeter & Koshland 1981)—does the model predict large-amplitude oscillations of the sort observed in the experiments.

Figure 8

Effect of zero-order ultrasensitivity in the phosphorylation–dephosphorylation kinetics of Msn2. The results have been obtained for K9=K10=K11=K12=(a) 0.1, (b) 0.5 and (c) 5, in the deterministic (left column) and the stochastic (right column) versions of the model. Moreover, Ω=75 for the stochastic case. Other parameters are the same as given in the legend to table 1.

Because phosphorylation–dephosphorylation cycles are also present in the part of the model devoted to the cAMP–PKA pathway, the question arises as to the effect of zero-order ultrasensitivity on oscillations of cAMP and, consequently, on Msn2. We thus tested the effect of increasing cycle by cycle the Michaelis constants that characterize the phosphorylation–dephosphorylation cycles acting on GEF, GAP, RAS, and PDE. The results indicate that for cAMP oscillations to occur, a sharp threshold is needed only in the cycle acting on RAS. Indeed, as long as this threshold exists, we may increase the values of parameters K1=K2, K3=K4 or K7=K8, one pair at a time, without affecting much the oscillations of cAMP. The role of the threshold in the cycle acting on the interconversion between RGDP and RGTP is illustrated in figure 9. When holding the other parameters at their basal values listed in the legend to table 1, we see that cAMP oscillations (a) and the associated oscillations of Msn2 (b) disappear when K5=K6 rise above a value close to 0.006×RASt (in μM). We further observe that oscillations appear through a supercritical Hopf bifurcation as the parameters K5=K6 decrease below the value 0.0022×RASt (in μM). This creates a small-amplitude stable limit cycle, which soon turns around and becomes unstable in a region where a large, stable limit cycle exists. Thus a stable steady state coexists with a stable limit cycle over a sizeable range of K5=K6 values, whereas near the bifurcation point two stable limit cycles of large or small amplitude coexist over a narrow range of K5=K6 values, a situation known as birhythmicity (Decroly & Goldbeter 1982).

Figure 9

Role of zero-order ultrasensitivity in the origin of cAMP oscillations. The bifurcation diagram shows the envelope of oscillations for (a) cAMP and (b) nuclear Msn2, both expressed in μM, as a function of K5=K6 (divided by RASt) as well as the stable steady state (middle branch) in the deterministic model. Other parameter values are the same as given in the legend to table 1. Dots representing data obtained by numerical integration of equations (2.1a)–(2.1g) and (2.3a)–(2.3d) correspond to stable steady-state or periodic behaviour. Solid and dashed lines have been obtained with Auto (Doedel 1981) and are associated, respectively, with stable and unstable periodic behaviours. The Hopf bifurcation is located at K5=K6=0.00215×RASt. Two limit points have been identified at K5=K6=0.00205×RASt, and at K5=K6=0.0057×RASt. Between these two values we observe coexistence between a stable steady state and a stable limit cycle.

The effect of zero-order ultrasensitivity in the RAS interconversion cycle can be assessed by determining the time evolution as a function of increasing values of K5=K6, in both the deterministic and stochastic versions of the model. A most conspicuous property brought to light by numerical simulations (figure 10) is that stochastic time series show undamped oscillations for cAMP and Msn2 in both figure 10a and 10b, even though the values of K5=K6 correspond, respectively, to sustained and rapidly damped oscillations in the deterministic case. This result is obtained even further away from the oscillatory domain, for example for K5=K6=0.1×RASt. It is only when zero-order ultrasensitivity in the RAS cycle disappears (figure 10c) that the stochastic time series becomes highly noisy and ceases to display high-amplitude oscillations.

Figure 10

Large-amplitude stochastic oscillations can occur outside the domain of sustained deterministic oscillations. The figure compares the time evolution of cAMP and nuclear Msn2 in the deterministic (left column) and stochastic (right column) versions of the model for increasing values of K5=K6: (a) 0.005×RASt, (b) 0.01×RASt and (c) 1.0×RASt. Stochastic simulations have been performed for Ω=75. Other parameter values are the same as given in the legend to table 1. cAMP and nuclear Msn2 are expressed in μM and in numbers of molecules in the deterministic and stochastic cases, respectively. In the presence of molecular noise, undamped large-amplitude oscillations are observed over a large range of K5=K6 values corresponding to rapidly damped oscillations in the deterministic case (see text).

The question arises as to whether the occurrence of large-amplitude stochastic oscillations outside the domain of deterministic oscillations is related to the existence of a region of hard excitation in which a stable steady state coexists with a stable, large-amplitude limit cycle in the diagram of figure 9. The mechanism of cAMP oscillations contains two negative feedback loops, both mediated through the catalytic subunit C of PKA (see figure 1): the first feedback is based on the activation of GAP, while the second is based on the activation of PDE. Each of these feedback loops is capable of producing sustained oscillations of cAMP (Garmendia-Torres et al. 2007). When we suppress the feedback on PDE and retain only the negative feedback exerted via GAP, then the bifurcation diagram as a function of K5=K6 loses its domain of hard excitation. For the parameter values considered in figure 9, the oscillations progressively decrease in amplitude and finally disappear as K5=K6 rise above the critical value close to 0.11×RASt (see figure S1 in the electronic supplementary material). This shows first that oscillations do not always require strong zero-order ultrasensitivity in RAS modification. Second, when stochastic simulations are performed to the left and to the right of the bifurcation point, e.g. for K5=K6=0.05×RASt or 0.2×RASt (see figure S2 in the electronic supplementary material), we do not observe the occurrence of large-amplitude stochastic oscillations outside the deterministic oscillatory domain, as was the case in figure 9.

7. Discussion

Recent experimental studies have shown that stress in yeast induces nucleocytoplasmic oscillations of the transcription factor Msn2. Corroborated by the study of a deterministic computational model, subsequent experiments suggested that the oscillations of Msn2, of a period of the order of minutes, originate from stress-induced oscillations of cAMP. The latter are accompanied by a periodic variation in PKA activity, which, in turn, controls the oscillatory changes in subcellular localization of Msn2. Fluorescence measurements of Msn2–GFP in individual yeast cells show that the oscillations in Msn2 shuttling are highly irregular. If the deterministic model can account for the oscillatory dynamics of Msn2 shuttling, it is however unable to reproduce the irregular oscillatory behaviour seen in the experiments. We therefore developed in this work a stochastic approach to determine the role of molecular noise in the origin of irregular Msn2 oscillations within individual yeast cells.

We first transformed the deterministic model described by equations (2.1a)–(2.3d) into a non-developed stochastic version containing 22 reaction steps (table 1). The stochastic oscillations predicted by this model (figure 3) were compared with the experimental observations on oscillatory Msn2 shuttling (figure 4). We determined the effect of a decrease in the number of molecules on the robustness of stochastic oscillations in Msn2 and cAMP (figures 5 and 6). Finally, we showed that the robustness of oscillations with respect to molecular noise depends on the distance from the bifurcation points (figure 7) determined in the deterministic model (figure 2). The stochastic study indicates that large-amplitude oscillations of cAMP are reflected by large-amplitude oscillatory changes in Msn2 subcellular localization only in the middle of the range of stress intensity corresponding to limit cycle behaviour.

Owing to the uncertainty on parameter values, it is difficult at this stage to make specific predictions as to the number of molecules needed to ensure robust oscillations. However, the results of figure 5 confirm the trend to be expected when lowering the number of molecules involved in the oscillatory mechanism. The results further show that coherent oscillations are progressively overcome by noise when the total number of Msn2 molecules decreases from 75 to 25. Moreover, the results of figure 8 suggest that at least moderate zero-order ultrasensitivity in Msn2 reversible phosphorylation is likely to be present in vivo to account for the sizeable periodic variations observed for Msn2 shuttling in intact yeast.

In regard to the role of zero-order ultrasensitivity in the cAMP–PKA pathway, a relatively steep threshold in the modification cycle involving RGDP and RGTP appears to be needed for sustained oscillations. This is not the case for the other proteins such as GEF, GAP, and PDE. Steep ultrasensitivity in the RAS switch appears to be favoured by the large levels of RAS observed experimentally (Garmendia-Torres et al. 2007). In determining the role of zero-order ultrasensitivity, we established the bifurcation diagram as a function of the Michaelis constants K5 and K6, which are associated with the RAS interconversion cycle. We found that a stable steady state coexists with a stable limit cycle at the extremity of the oscillatory domain, a phenomenon known as hard excitation. Unexpectedly stochastic simulations showed in these conditions the occurrence of high-amplitude oscillations well outside the domain of sustained oscillations predicted by the deterministic model. This behaviour, which represents an instance of noise-induced oscillations (Horsthemke & Lefever 1983), occurs close to a domain of hard excitation, and is not seen in the absence of such a domain (see figures S1 and S2 in the electronic supplementary material).

In coupling Msn2 shuttling to PKA oscillations, we considered for simplicity a single phosphorylation of the protein by PKA. We know that four PKA sites control the NLS, and the NES also contains a PKA phosphorylation site. It appears that more than two sites of the NLS need to be dephosphorylated to activate it (C. Garmendia-Torres and M. Jacquet 2007, unpublished observations). Multiple phosphorylations may lead to more complex regulation, as is well illustrated in the case of the yeast transcription factor Pho4 (Springer et al. 2003). This transcriptional activator displays different degrees of phosphorylation as a function of the amount of inorganic phosphate in the medium. At a high phosphate level, Pho4 is fully phosphorylated and remains cytoplasmic. Under phosphate limitation there is only one specific phosphorylation; Pho4 then goes to the nucleus and activates a set of genes. In the absence of phosphate in the medium, no phosphorylation occurs; Pho4 then interacts with another effector and regulates another set of genes. It is likely that the multiple sites of phosphorylation of Msn2 by PKA could generate also a more complex behaviour than the one presented here.

The reason why we consider a single phosphorylation is that, as stated above in §2.1, this suffices to couple Msn2 shuttling to cAMP oscillations and PKA. Moreover, the kinetics and role of multiple phosphorylation of Msn2 are far from being known in as much detail as in the case of Pho4 (Jeffery et al. 2001). Incorporating multiple phosphorylation of Msn2 into the deterministic and stochastic models for oscillatory shuttling of this transcription factor represents a natural extension that awaits further experimental information.

Stochastic simulations of the oscillatory nucleocytoplasmic shuttling of Msn2 coupled to oscillations in the cAMP–PKA pathway yield good agreement with the experimentally observed oscillations of Msn2. This supports the view that the irregular nature of the oscillations in nucleocytoplasmic shuttling of Msn2 originates from the effect of molecular noise at the level of single yeast cells. Molecular noise is associated with the relatively small numbers of molecules involved in the mechanism of oscillations underlying the oscillatory shuttling of Msn2.

Much work has been devoted to the characterization of the effects of noise in gene expression, in both bacteria (Elowitz et al. 2002; Ozbudak et al. 2002; Süel et al. 2007) and eukaryotes (Blake et al. 2003; Raser & O'Shea 2004, 2005). In complement to experimental studies in single cells, the role of fluctuations in gene transcription and in genetic networks has been investigated theoretically (Arkin et al. 1998; Swain et al. 2002; Swain 2004; Scott et al. 2006). In regard to oscillatory behaviour, the effect of noise has been considered in theoretical studies of genetic oscillators such as those underlying circadian rhythms (Barkai & Leibler 2000; Gonze et al. 2002a, 2004; Vilar et al. 2002; Leloup et al. 2006) or the Repressilator, a synthetic gene circuit displaying oscillatory transcriptional activity (Elowitz & Leibler 2000; Scott et al. 2006).

Here we considered the effect of noise in a model for a cellular oscillatory process based on metabolic rather than genetic regulation. The comparison between the results of deterministic and stochastic simulations of periodic Msn2 shuttling controlled by oscillations in the cAMP–PKA pathway indicates that the deterministic approach provides us with a qualitative picture that encompasses well the most conspicuous properties of this new oscillatory phenomenon, besides its irregular nature due to stochastic fluctuations. Beyond the difference in the nature of the underlying regulatory process, this conclusion holds with that reached (Gonze et al. 2002a,b, 2004; Leloup et al. 2006) when comparing stochastic and deterministic models for circadian rhythms based on genetic regulation and post-transcriptional modification.

Acknowledgments

We wish to thank the referees for their fruitful suggestions. The work was supported by grant no. 3.4636.04 from the Fonds de la Recherche Scientifique Médicale (F.R.S.M., Belgium), by the European Union through the Network of Excellence BioSim, Contract no. LSHB-CT-2004-005137, and by the Belgian Programme on Interuniversity Attraction Poles, initiated by the Belgian Federal Science Policy Office, project P6/25 (BioMaGNet). This work was initiated while the authors took part in the workshop ‘Biological switches and clocks’ organized at KITP, University of California at Santa Barbara, in July 2007.

Footnotes

References

View Abstract