## 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.

## 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*. 2002*a*).

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).

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):(2.1a)(2.1b)(2.1c)(2.1d)(2.1e)(2.1f)(2.1g)with(2.2a)

(2.2b)

(2.2c)

In the above equations, GEF_{t}, GAP_{t}, RAS_{t}, CYCL_{t}, PDE_{t} and PKA_{t} 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.1*b*) and (2.1*e*) by assuming that the catalytic subunit C of PKA catalyzes the phosphorylation of the inactive GAP protein, GAP_{i}, into the active form GAP_{a}, 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:(2.3a)(2.3b)(2.3c)(2.3d)where *k*_{t1}, *k*_{t2}, *k*_{t3} and *k*_{t4} denote first-order apparent rate constants for transport of the various forms of Msn2 into and out of the nucleus. *V*_{9} and *V*_{11} 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.1*g*), (2.2*a*)–(2.2*c*) 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.3*b*) and (2.3*c*), 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.1*f*). 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 *V*_{10} and *V*_{12} are multiplied by parameter *Str*.

The total concentrations of Msn2 in the cytosol (*M*_{cyto}) and in the nucleus (*M*_{nucl}) are given by(2.4a)

(2.4b)

The conservation equation MSN2_{t}=*M*_{cyto}+*M*_{nucl} 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.3*a*)–(2.3*d*).

Equations (2.1*a*)–(2.4*b*) 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 2*a*), which leads to the predominantly cytoplasmic localization of Msn2 (figure 2*b*). Beyond a critical value of stress intensity, oscillations of cAMP occur spontaneously (figure 2*a*). 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 3*a*). Finally, above a second, higher critical value of stress intensity, cAMP reaches a stable, low steady-state level (figure 2*a*), while Msn2 ceases to oscillate and becomes predominantly located in the nucleus (figure 2*b*).

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 2*b*).

## 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*. 2002*b*, 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 *p*_{i} of reaction step *i* is then calculated from the reaction propensity *w*_{i} according to the expression .

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 3*b* 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 2*a* 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 3*a*). 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 4*a* are typical oscillations of nuclear (and, in antiphase, cytoplasmic) Msn2 observed in the experiments. Figure 4*b* 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 4*a* and 4*b* 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 4*c* 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.

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 (*b*–*d*) 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 5*a*–*d* (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, MSN_{t}, which is taken as equal to 1 μM (see legend to table 1). The comparison of the curves in figure 5*c*,*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.

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 5*d* (also established for *Ω*=25). The oscillations in cAMP (figure 6*a*) are indeed as robust as those shown in figure 5*b* (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.

## 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 2*a*. 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 7*a*). 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 7*b*, established for *Str*=0.75). When the system is located well inside the oscillatory domain (*Str*=1), periodic behaviour is clearly discernible (figure 7*c*) and the phase plane trajectory takes the form of a high-amplitude, noisy limit cycle surrounding the deterministic limit cycle trajectory (white curve).

For a higher value of stress intensity (*Str*=1.3) the deterministic bifurcation diagram of figure 2*a* 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 7*c* and figure 7*d*). 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 7*e*).

## 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 7*c*,*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 *K*_{9}–*K*_{12} 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 8*a*)—which conditions favour ultrasensitivity (Goldbeter & Koshland 1981)—does the model predict large-amplitude oscillations of the sort observed in the experiments.

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 *K*_{1}=*K*_{2}, *K*_{3}=*K*_{4} or *K*_{7}=*K*_{8}, 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 *K*_{5}=*K*_{6} rise above a value close to 0.006×RAS_{t} (in μM). We further observe that oscillations appear through a supercritical Hopf bifurcation as the parameters *K*_{5}=*K*_{6} decrease below the value 0.0022×RAS_{t} (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 *K*_{5}=*K*_{6} values, whereas near the bifurcation point two stable limit cycles of large or small amplitude coexist over a narrow range of *K*_{5}=*K*_{6} values, a situation known as birhythmicity (Decroly & Goldbeter 1982).

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 *K*_{5}=*K*_{6}, 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 10*a* and 10*b*, even though the values of *K*_{5}=*K*_{6} 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 *K*_{5}=*K*_{6}=0.1×RAS_{t}. It is only when zero-order ultrasensitivity in the RAS cycle disappears (figure 10*c*) that the stochastic time series becomes highly noisy and ceases to display high-amplitude oscillations.

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 *K*_{5}=*K*_{6} loses its domain of hard excitation. For the parameter values considered in figure 9, the oscillations progressively decrease in amplitude and finally disappear as *K*_{5}=*K*_{6} rise above the critical value close to 0.11×RAS_{t} (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 *K*_{5}=*K*_{6}=0.05×RAS_{t} or 0.2×RAS_{t} (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.1*a*)–(2.3*d*) 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 *K*_{5} and *K*_{6}, 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*. 2002*a*, 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*. 2002*a*,*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

Electronic supplementary material is available at http://dx.doi.org/10.1098/rsif.2008.0141.focus or via http://journals.royalsociety.org.

One contribution of 10 to a Theme Supplement ‘Biological switches and clocks’.

- Received April 10, 2008.
- Accepted April 25, 2008.

- © 2008 The Royal Society