Multiscale Modeling of Bistability in the Yeast Polarity Circuit

TL;DR

This paper uses multiscale modeling to study how yeast cells break symmetry and establish polarity, revealing mechanisms like bistability and rapid dismantling of polarity sites.

Contribution

The novelty lies in combining ODE and RDE models across scales to explain symmetry breaking and bistability in yeast polarity.

Findings

A model with positive and negative feedback on different time scales explains bistability and rapid dismantling.

Coarse-grained RDE models capture spatial profiles of polarity factors effectively.

Particle-based models reveal stochastic fluctuations in biochemical reactions.

Abstract

Cell polarity refers to the asymmetric distribution of proteins and other molecules along a specified axis within a cell. Polarity establishment is the first step in many cellular processes. For example, directed growth or migration requires the formation of a cell front and back. In many cases, polarity occurs in the absence of spatial cues. That is, the cell undergoes symmetry breaking. Understanding the molecular mechanisms that allow cells to break symmetry and polarize requires computational models that span multiple spatial and temporal scales. Here, we apply a multiscale modeling approach to examine the polarity circuit of yeast. In addition to symmetry breaking, experiments revealed two key features of the yeast polarity circuit: bistability and rapid dismantling of the polarity site following a loss of signal. We used modeling based on ordinary differential equations (ODEs) to investigate mechanisms that generate these behaviors. Our analysis revealed that a model involving positive and negative feedback acting on different time scales captured both features. We then extend our ODE model into a coarse-grained reaction–diffusion equation (RDE) model to capture the spatial profiles of polarity factors. After establishing that the coarse-grained RDE model qualitatively captures key features of the polarity circuit, we expand it to more accurately capture the biochemical reactions involved in the system. We convert the expanded model to a particle-based model that resolves individual molecules and captures fluctuations that arise from the stochastic nature of biochemical reactions. Our models assume that negative regulation results from negative feedback. However, experimental observations do not rule out the possibility that negative regulation occurs through an incoherent feedforward loop. Therefore, we conclude by using our RDE model to suggest how negative feedback might be distinguished from incoherent feedforward regulation.

Figures

Full text

1. Introduction

Cells are constantly responding to changes in their environment [1,2,3]. Often, the proper response to an environmental cue requires the generation or reorganization of internal structures. For example, yeast cells undergo directed growth in response to gradients of mating pheromone [4,5,6]. This chemotropic growth requires the formation of an internal structure referred to as the polarity site [7]. The first step in the formation of the polarity site is the establishment of a localized region of high Cdc42 concentration at the cell periphery [8,9,10,11]. Cdc42 is a member of the Rho family of GTPases and functions as a molecular switch, existing in either an inactive or active state [12,13,14]. In its active state, Cdc42 is located in the cell membrane and recruits actin cables along which vesicles are trafficked. The insertion of these vesicles into the cell membrane generates growth in the direction of the polarity site. The biochemical reaction network that generates clusters of Cdc42 activity is referred to as the polarity circuit (Figure 1A). Polarity establishment occurs via a positive feedback loop in which active Cdc42 recruits its activator, the GEF Cdc24, in complex with the scaffold Bem1 from the cytosol. In turn, Cdc24 activates neighboring Cdc42 molecules. This positive feedback loop, coupled with the different diffusion rates associated with membrane-bound active Cdc42 (slow diffusion) and cytosolic inactive Cdc42 (rapid diffusion), is sufficient to generate a polarity site through a Turing mechanism [15]. The ability of the yeast polarity circuit to undergo Turing instability was first demonstrated through mathematical modeling [16]. Since then, mathematical modeling has been used to explain many properties of the polarity circuit.

The kinetics of biochemical networks can be modeled using rate equations that take the form of ordinary differential equations (ODEs) [17,18,19,20]. ODEs represent a coarse-grained description of the system, are computationally efficient to simulate, and often lend themselves to systematic analysis. However, they do not allow for spatial variations in concentrations, which are critical for understanding cell polarity. Therefore, the most common way to model the polarity circuit is using reaction–diffusion equations (RDEs) that take the form of partial differential equations [21,22,23,24,25]. RDEs capture both the biochemical reactions involved in the polarity circuit and the thermal diffusion of molecules. The drawbacks of using RDEs are that they can be computationally expensive to simulate, and they do not capture fluctuations arising from the stochastic nature of biochemical reactions. The dynamics of biochemical networks are most accurately captured by particle-based models [26,27,28]. In these models, every molecule in the system is considered. Molecules undergo thermal diffusion and can react if they are sufficiently close to a reaction partner. These simulations account for stochastic effects due to diffusion and the inherently random nature of chemical reactions. The drawbacks of particle-based models are that they are computationally expensive to simulate and do not lend themselves to theoretical analysis.

We apply each modeling approach to investigate key properties of the yeast polarity circuit across spatial scales of increasing resolution. Yeast cells undergo polarized growth in two different contexts: budding and mating [29,30,31,32]. Polarization during mating requires the same positive feedback circuit that is crucial for polarization during budding. In contrast to budding, polarity establishment during mating requires an external stimulus, mating pheromone, that can easily be manipulated under laboratory conditions. Using an ODE model of the pheromone-induced polarity circuit, we demonstrate that adding negative regulation to the core polarity circuit is sufficient to capture several dynamic properties of yeast polarity establishment, including bistability and the rapid disassembly of the polarity site following the removal of pheromone. We then converted the ODE model to an RDE model by adding terms that account for diffusion. After verifying that the RDE model performs as expected, we extended the model to more faithfully capture the biochemical reactions involved in the polarity circuit. Finally, we recast our detailed RDE model as a particle-based model to demonstrate it still performs as expected when stochastic effects due to the random nature of biochemical reactions are considered. All our models assume that negative regulation of the polarity site occurs through negative feedback. However, it is also possible that an incoherent feedforward loop underlies the rapid dismantling of the polarity site. Therefore, we end by performing an analysis of the RDE model to demonstrate how, in theory, negative feedback can be distinguished from incoherent feedforward regulation.

2. Methods

2.1. Experimental Datasets

Files containing the experimental data referenced in this manuscript and a document containing an index for the datasets and an explanation of the image analysis methods used to quantify polarity are available on dryad (https://doi.org/10.5061/dryad.4f4qrfjn0).

2.2. Dose Response Curves Derived from Experimental Data

The polarity status as pheromone concentration decreases from 50 nM to 0 nM was derived from the dataset indexed as Experiment 7 (https://doi.org/10.5061/dryad.4f4qrfjn0). The average deviation of uniformity (DfU) for the three time points at each pheromone concentration was calculated. The averages and standard error were calculated from five independent experiments. The polarity status as pheromone concentration increases from 0 nM to 50 nM was derived from datasets indexed as Experiments 1–5 (https://doi.org/10.5061/dryad.4f4qrfjn0) as follows. The average DfU for 0 nM and 0 to 4 nM pheromone concentrations were calculated for time points across the 55 min and 70 min time course, respectively, because polarization was not achieved (DfU < 0.35). The average DfU for 6–50 nM was calculated from all time points in the 120 min time course after a polarized status was achieved (DfU > 0.4). The average DfU for 5 nM pheromone was calculated from the last seven time points in the time course. The averages and standard error were calculated from three independent experiments.

2.3. Reaction–Diffusion Equation Simulations

Both the simplified and detailed RDE models were solved in 1D using periodic boundary conditions. To generate the time series presented in the figures, we numerically solved the equations using a forward Euler method with a domain length of $[eqn]$ . The equations were solved using a spatial step size of $[eqn]$ . The time step was initially set to $[eqn]$ s and reduced until the solution no longer depended on $[eqn]$ . The bifurcation diagrams for RDE models were generated by simulating the system starting from (1) homogeneous initial conditions to which a small amount of noise was added and (2) a localized initial spike of Cdc42 activity. Following the simulations, the maximum value of active Cdc42 concentration was recorded for each initial condition and used to determine boundaries between the monostable (low active Cdc42 for both initial conditions), bistable (low active Cdc42 for the uniform initial condition and high for the spike initial condition), and Turing unstable (high active Cdc42 for both initial conditions) regimes. To create the bifurcation diagrams, we solved the equations using the MATLAB (R2023b) solver pdepe with no-flux boundary conditions and a domain size of $[eqn]$ . At a steady state, a solution with periodic boundary conditions can be constructed by reflecting the no-flux solution around the boundary where polarity occurred. We selected several parameter sets to verify that the two simulations methods produced identical steady state solutions.

Parameters for the simplified model are listed in Table S1, and those for the detailed model are listed in Table S2. Most reaction rates for the detailed RDE model were scaled from those for the 2D RDE model in Pablo et al. [28]. All code is available at https://github.com/guanhighyun/yeast-polarity-negative-feedback (accessed on Aug 13 2024).

2.4. Particle-Based Simulations

Particle-based simulations were performed using Smoldyn (v2.67) on a Linux-based computing system (Longleaf cluster at UNC, Chapel Hill, NC, USA) [27]. The code is available at https://github.com/guanhighyun/yeast-polarity-negative-feedback (accessed on 13 August 2024). Molecules were regarded as point particles and conducted Brownian motion according to the Euler–Maruyama method. We used a 2D domain with length = $[eqn]$ and periodic boundary conditions for all simulations. Rates for first-order and second-order inter-particle reactions are listed in Table S3. All of the parameters referred to a previous publication [28] except for GAP-associated reactions ( $[eqn]$ , $[eqn]$ , and $[eqn]$ ). Polarity strength was measured using Ripley’s K-function ( $[eqn]$ ) [28,33]. $[eqn]$ when particles are distributed uniformly. A high $[eqn]$ value ( $[eqn]$ ) signifies strong polarization.

3. Results

3.1. Experimental Observations

Data demonstrating that the yeast polarity circuit is bistable with respect to pheromone concentration was presented by Vered [34]. These observations were made using a microfluidics chamber to precisely control the temporal profile of pheromone and live-cell imaging of fluorescent reporters to monitor the polarity process in single living cells. Details of the microfluidic device used in these experiments have been published elsewhere [35,36,37]. Here, we summarize observations from three sets of experiments that provided insights into the circuit that were foundational for our model development. (The associated datasets along with an overview of the experiments, description of the experimental methods, and image analysis techniques used to quantify polarity are published on DRYAD: https://doi.org/10.5061/dryad.4f4qrfjn0).

One set of experiments monitored polarity establishment upon exposure to pheromone concentrations of 4, 5, 6, 10, and 50 nM. Analysis of these datasets showed that 6 nM pheromone is the minimum concentration required to establish and maintain a stable polarity patch (See Figure 3.1A–C of [34]). Another notable observation was that at 5 nM pheromone, some cells failed to polarize, but others established a weak and unstable polarity patch that disassembled and reassembled repeatedly throughout the experiment (See Figure 3.1B of [34] and Movie S1).

In another experiment, polarity was established by exposing cells to 50 nM pheromone for 2 h, dropping the concentration to 10 nM (where polarity is stable), and then decreasing the concentration to 0 in 1 nM increments (See Figure 3.2A of [34]). Under this regimen, most cells continued to hold polarity at concentrations below the minimum pheromone concentration required for naïve cells to establish polarity, with some cells maintaining polarity even after the pheromone was completely removed.

Dose–response curves from these two datasets are shown in Figure 1B (see Methods for details for how these curves were constructed and https://doi.org/10.5061/dryad.4f4qrfjn0 for a discussion of using the deviation from uniformity (DfU) as a measure of polarity). The curves demonstrate that polarization as a function of pheromone concentration depends on the history of the system. There is a range of pheromone concentrations for which two steady states co-exist (i.e., the system is bistable). How the system is prepared determines which of these two steady states is experimentally observed (i.e., the system exhibits hysteresis). If we start from low levels of pheromone and increase the concentration, the system will be in the unpolarized state of the bistable regime (Figure 1B, black curve). Alternatively, if we start from high levels of pheromone and gradually reduce the concentration, the system will be in the polarized state of the bistable regime (Figure 1B, red curve). Another way to state this hysteretic effect is that the minimum pheromone concentration required to establish polarity is greater than the minimum pheromone concentration required to maintain polarity. Such hysteretic behavior is a hallmark of bistability.

In another series of experiments, cells were exposed to a high pheromone concentration (50 nM) to establish polarity, and then the pheromone level was abruptly reduced to 0, 4, 5, 6, or 10 nM. Polarity was lost within 5 min when the pheromone was reduced to 5 nM or less (See Figure 3.3A,B of [34]). The same polarity loss occurs for cells exposed to intermediate pheromone concentrations (10 nM) before the pheromone level is abruptly dropped to 0 (see Figure 3.5B of [34]). By contrast, cells polarized by exposure to intermediate pheromone levels (10 nM) can maintain the polarity patch if the pheromone is not completely removed but abruptly dropped to 5 nM (see Figure 3.5B of [34]).

3.2. Mechanistic Considerations for Model Development

We reasoned two mechanisms could account for the rapid loss of polarity. In the first scenario, inhibition in the form of either negative feedback (Figure 2A) or incoherent feedforward regulation (Figure 2B) dismantles the polarity site following the removal of the pheromone. In the second mechanism, the slow reinforcement of the polarity site acts to stabilize the site over time. In this scenario, polarity is lost when insufficient time has passed for reinforcement to occur before the pheromone is removed. The reinforcement mechanism is motivated by the observation that actin cables are required for long-term maintenance of the polarity site but not its initial formation [38]. However, our experimental observation that a drop from 10 to 5 nM pheromone levels does not cause a loss of polarity argues against this scenario. Reinforcement would be the same or diminished compared with the drop from 50 to 5 nM, where polarity is destroyed. Therefore, we focus on the negative regulation mechanism. To test if negative regulation is sufficient to explain our experimental results, we developed a series of models of increasing spatial resolution. We focus on negative feedback as the regulatory mechanism responsible for dismantling the polarity site but note that currently available data do not allow us to rule out incoherent feedforward regulation. We end our analysis by suggesting an experiment that might distinguish the two negative regulation mechanisms.

3.3. Ordinary Differential Equation Model

The observed rapid disassembly of the polarity patch following the removal of the pheromone suggests the presence of pheromone-induced negative regulation that remains active after the pheromone is removed. In particular, it implies that following pheromone removal, this negative regulation is strong enough to disrupt the positive feedback loop that establishes polarity. By contrast, when the pheromone level is reduced slowly, our results suggest that positive and negative regulation remain balanced so that polarity is maintained. As an initial test of these ideas, we developed an ODE model that consists of a positive feedback loop that generates bistability and a slow negative feedback loop. In the model, Cdc42 exists in an active GTP bound state ( $[eqn]$ ) or an inactive GDP bound state ( $[eqn]$ ). We assume high levels of $[eqn]$ represent a polarized state. We also assume that the total amount of Cdc42 ( $[eqn]$ ) is constant, so that $[eqn]$ . The equation that governs the time evolution of $[eqn]$ takes the following form:

$$[eqn]$$

The first term on the right-hand side in the above equation models the conversion of inactive Cdc42 to its active form. The parameter $[eqn]$ is the basal activation rate. The term $[eqn]$ models the positive feedback loop in which active Cdc42 promotes the activation of inactive Cdc42. Biologically, this term represents the recruitment of the GEF Cdc24 by active Cdc42. We assume that the effect of pheromone is to increase this rate: $[eqn]$ , where $[eqn]$ is the basal (pheromone independent) rate and $[eqn]$ represents the increase in this rate in the presence of the pheromone. This would occur, for example, if pheromone-induced MAPK activity increased the rate at which the GEF for Cdc42 is recruited to the membrane. The second term in the above equation represents the rate at which active Cdc42 becomes inactive. The parameter $[eqn]$ is the basal rate of inactivation, and the term $[eqn]$ models the effect of an (as yet) unidentified negative regulator ( $[eqn]$ ). We assume the equation that governs the time evolution of the negative regulator is given by the following:

$$[eqn]$$

The first term on the right-hand side of the above equation represents the rate at which $[eqn]$ increases in response to active Cdc42 ( $[eqn]$ ), and the second term represents the rate at which the negative regulator ( $[eqn]$ ) is lost. Our model represents a well-stirred version of the model considered by Chiou et al. [39] to study the competition and coexistence of polarity sites with the additional feature of a negative feedback loop.

For bistability, Equations (1) and (2) must admit multiple steady-state solutions. A steady state refers to the condition in which $[eqn]$ and $[eqn]$ are constant in time. That is, $[eqn]$ . Therefore, to find the steady state values of $[eqn]$ and $[eqn]$ , we solve the following set of equations: $[eqn]$ and $[eqn]$ . The two equations define curves in the $[eqn]$ phase plane referred to as nullclines (Figure 3A). The steady states of the system occur where the nullclines intersect. When $[eqn]$ (no pheromone), there are three steady states shown in Figure 3A as intersection of the blue line ( $[eqn]$ -nullcline) and black curve ( $[eqn]$ -nullcline). The intersection points at low and high values of $[eqn]$ represent the unpolarized and polarized states of the system, respectively. The steady state at the intermediate value of $[eqn]$ is unstable and cannot be observed experimentally. Therefore, in the absence of pheromone, the system is bistable. When $[eqn]$ is increased beyond a critical point, the lower two steady states are lost, and the system transitions to a monostable polarized state, shown in Figure 3A as the intersection of the blue line ( $[eqn]$ nullcline) and the solid red curve ( $[eqn]$ nullcline labeled with “ $[eqn]$ ”).

To investigate the system’s temporal response to pheromone, we start with $[eqn]$ and $[eqn]$ and X at their steady-state values in the unpolarized state. We then turn on pheromone by setting $[eqn]$ and numerically solve Equations (1) and (2) (Figure 3B). The concentration of $[eqn]$ (active Cdc42) rapidly increases as the system polarizes. This rapid rise in $[eqn]$ is followed by a gradual decline as the level of $[eqn]$ (negative regulator) increases toward its steady-state value. To gain further insight into the behavior of the system, we can plot $[eqn]$ versus $[eqn]$ in the $[eqn]$ plane (Figure 3A, dashed red curve labeled “signal on”). As can be seen from this graph, $[eqn]$ rapidly approaches the $[eqn]$ nullcline for $[eqn]$ and then follows this curve to the fixed point as $[eqn]$ slowly increases. We next use Equations (1) and (2) to simulate the removal of pheromone by starting the model from the steady state for $[eqn]$ and then setting $[eqn]$ (Figure 3C). In this scenario, $[eqn]$ rapidly decreases, indicating a loss of polarity, while the value of $[eqn]$ decreases over a longer time scale. Again, plotting $[eqn]$ versus $[eqn]$ in the $[eqn]$ plane, we see that $[eqn]$ rapidly approaches the lower branch of the $[eqn]$ nullcline, before $[eqn]$ significantly changes (Figure 3A, dashed red curve labeled “signal off”). The system then follows this nullcline to the unpolarized steady state as $[eqn]$ (negative regulator) slowly decreases. Finally, we simulated a drop in pheromone to a level of $[eqn]$ from two different starting concentrations (Figure 3D). When the initial value of $[eqn]$ was taken to be $[eqn]$ the system lost polarity when $[eqn]$ was dropped to $[eqn]$ (Figure 3D, red curve). However, when the initial value of $[eqn]$ was used, polarity was maintained following the drop to $[eqn]$ (Figure 3D, green curve).

These results provide evidence in support of a bistable system coupled with slow negative feedback. However, the ODE model does not consider spatial aspects of polarity establishment, and therefore, it is not clear if the results still apply when spatiotemporal models are considered.

3.4. Reaction–Diffusion Equation Models

Next, we updated our model to include spatial variations in the concentrations of the molecular species in the polarity circuit. To start our analysis, we considered a model without negative feedback:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is localized to the membrane and $[eqn]$ is located in the cytosol. The last terms in Equations (3) and (4) account for the effects of diffusion. Importantly, diffusion in the membrane is substantially slower than diffusion in the cytosol (i.e., $[eqn]$ ).

To begin our analysis, we constructed a two-parameter bifurcation diagram using the rate constants $[eqn]$ (positive feedback) and $[eqn]$ (deactivation) to characterize the steady-state behavior of Equations (3) and (4) (Figure 4A). This diagram consists of three regions. In the bistable region (white area), the spatially homogenous and polarized steady states are stable to small perturbations. In the Turing unstable region (red area), the spatially homogenous solution is unstable, and the system always evolves to a polarized state. In the monostable region (black area), the system does not polarize, and the only steady-state solution is spatially homogenous. As we demonstrate below, when negative regulation is added to the model, we can interpret the behavior of the extended model in terms of movement through the parameter space.

Next, we update the model equations to include the effect of pheromone and a negative feedback loop (Figure 2A). The revised model equations are as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Again, we assume that the effect of the pheromone is to increase the strength of the positive feedback loop through the relation $[eqn]$ . We replace $[eqn]$ with the expression $[eqn]$ , where $[eqn]$ is the concentration of the negative regulator. In Equation (7), we assume that the production rate of $[eqn]$ depends on $[eqn]$ , forming a negative feedback loop. We also assume that $[eqn]$ exists in the cytosol and, therefore, diffuses at a rate governed by $[eqn]$ .

We are now able to use the bifurcation diagram to understand the behavior of the model with negative feedback. In the absence of pheromone, the system starts in the bistable regime in the unpolarized steady state (Figure 4A, position 1). The addition of pheromone increases $[eqn]$ and, therefore, moves the system into the Turing unstable regime, leading to polarization (Figure 4A, position 2). Next, there is a slow buildup of $[eqn]$ that increases $[eqn]$ (Figure 4A, position 3). If $[eqn]$ is reduced slowly, the activation rate, $[eqn]$ , and deactivation rate, $[eqn]$ , remain balanced, and the system stays polarized after pheromone is completely removed (Figure 4A, position 4). However, if $[eqn]$ is reduced rapidly to zero, the activation rate drops immediately, but the delayed response of the negative feedback loop allows $[eqn]$ to remain elevated long enough to destroy polarity, and the system returns to the spatially homogenous steady state (Figure 4A, position 1).

To support this qualitative explanation, we performed computer simulations of the RDE model given by Equations (5)–(7) in a 1D domain using periodic boundary conditions. Full details for the numerical methods are provided in the Methods section. We first constructed a one-parameter bifurcation diagram using $[eqn]$ as the bifurcation parameter (Figure 4B). Consistent with our experimental results (Figure 1B), the model is bistable when $[eqn]$ and undergoes a Turing bifurcation at $[eqn]$ . Next, we simulated the case in which the system is started in the Turing unstable regime ( $[eqn]$ ), and $[eqn]$ is slowly decreased in time starting at $[eqn]$ min. (Figure 4C). As can be seen, both the active Cdc42 concentration, $[eqn]$ , and $[eqn]$ concentration have time to respond to the slowly changing signal, and polarity is maintained. In the second simulation, we used the same initial condition ( $[eqn]$ ). However, this time, the signal was turned off at $[eqn]$ min (Figure 4D). In this scenario, the polarity site is dismantled after the signal is removed.

3.5. Detailed RDE Model

While the model above captures qualitative properties of the polarity circuit, it lacks many of the biochemical steps involved in the system. To test if a more accurate representation of the system would behave in a similar fashion, we started with the model and parameter values presented by Pablo et al. [28]. To this model, we added a negative feedback loop in which active Cdc42 increases the rate at which its GAP gets activated. A schematic diagram of the extended model is shown in Figure 5A. This system can be represented with the following set of reaction–diffusion equations:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

All the parameter values are given in Table S2.

We first verified that the system produces bistable behavior as the pheromone concentration (implicitly modeled in the positive feedback rate $[eqn]$ ) is increased (Figure 5B). We then generated a two-parameter bifurcation diagram for the system lacking negative feedback ( $[eqn]$ ) as performed above for the simpler model (Figure 5C). The qualitative features of the bifurcation diagram remained unchanged for the extended model. We verified that the full system would maintain polarity if the pheromone concentration was reduced slowly and lose polarity if pheromone was rapidly removed (Figure 5D,E). Overall, results from our detailed RDE model align closely with the outcomes of the simplified model, effectively encapsulating yeast polarity features with a more realistic representation of reaction networks.

3.6. Particle-Based Model

While RDE models capture many qualitative features of the yeast polarity circuit, one of their shortcomings is that they do not capture the time scale over which polarity occurs, taking longer to polarize than observed experimentally [28]. Therefore, we turned to particle-based simulations to investigate how molecular-level fluctuations affect the performance of the polarity circuit [28]. We found that these intrinsic fluctuations substantially decrease the time needed for the system to polarize. However, the previous model did not include negative feedback or operate in a bistable regime. Therefore, we converted the detailed RDE model given by Equations (8)–(15) into a particle-based model that we simulated using the software application Smoldyn (v2.67) [27]. We used a cluster score based on Ripley’s K-function as described in the Methods section as a measure of polarity. Similar to the detailed RDE model, the particle-based model of the polarity circuit also showed a bistable region as the Cdc42 deactivation rate $[eqn]$ was varied (Figure 6A). When GAP-mediated negative feedback was added to the polarity circuit, the system continued to display bistable behavior and hysteresis as the positive feedback strength was varied (Figure 6B). To vary the positive feedback strength in the particle-based model, we varied the parameter $[eqn]$ , which is the rate at which inactive Cdc42 is converted to active Cdc42 if a Cdc42-GTP-GEF heterodimer is within the reaction radius. Consistent with our experimental observations, suddenly reducing $[eqn]$ destroyed polarity (Figure 6C, Movie S3), while slowly reducing $[eqn]$ allowed polarity to be maintained (Figure 6D, Movie S4). Interestingly, for the case in which the pheromone level is suddenly dropped, and polarity is lost, the polarity site eventually reformed at another location (Figure 6C). We also observed this behavior in a subset of our experiments (Movie S2). The time scale for polarization by the particle-based model was consistent with our experimental observations (~10–70 min). Our particle-based simulations confirm that the results obtained from our RDE simulations remain valid in the presence of molecular-level noise and capture additional features of the experimental results.

3.7. Negative Feedback Versus Incoherent Feedforward Regulation

In the models presented above, we assumed that negative feedback (NFB) was responsible for dismantling the polarity site following removal of pheromone. However, negative regulation through an incoherent feedforward (IFF) loop might also underlie this behavior. An incoherent feedforward loop is formed when a pheromone leads to the activation of a negative regulator of Cdc42 independently of Cdc42 activity (Figure 2B). To convert the model to one regulated by an incoherent feedforward loop, the only change we made was to the equation for the negative regulator:

$$[eqn]$$

The rest of the equations remained unchanged, and all parameter values remained the same, except for those in the equation above (see Table S1). In Equation (16), the synthesis rate of the negative regulator $[eqn]$ depends directly on the signal strength $[eqn]$ , as opposed to the amount of active Cdc42 as in the case of NFB (see Equation (7)). Our previous results from NFB simulations hold true for IFF simulations: (1) the system undergoes a bifurcation from a bistable region to a single steady state as the signal $[eqn]$ is increased (Figure S1A); (2) polarity is maintained if the pheromone signal is gradually reduced (Figure S1B) but is lost if the signal is rapidly removed (Figure S1C). The similar behavior of the NFB and IFF models made us wonder if it was, in theory, possible to distinguish these two forms of negative regulation.

For the NFB model, the concentration of the negative regulator $[eqn]$ is dependent on the amount of active Cdc42, whereas, for the IFF model, $[eqn]$ is independent of Cdc42 and only depends on the pheromone concentration. For both models, in the absence of pheromone ( $[eqn]$ ), the system is in a bistable regime and starts from the unpolarized state. When the pheromone level is increased, both systems move into a regime where only the polarized state is stable. If the system is started in this regime and the pheromone level is slowly reduced so that polarity is maintained, the amount of $[eqn]$ will be different in the two models. For NFB regulation, $[eqn]$ will remain elevated because its level depends on active Cdc42. In contrast, for IFF regulation, $[eqn]$ will decrease because, in this model, the level of $[eqn]$ depends only on the signal strength and not the amount of active Cdc42. This difference leads to a potential experimental method for distinguishing the two models.

The experiment is to start with no pheromone ( $[eqn]$ ) but using two different initial conditions. In the first case, the system starts from an unpolarized state. In the second case, the system is started from the polarized state. Experimentally, this would require exposing cells to high pheromone concentrations, then slowly reducing the concentration to zero. Both initial conditions are exposed to a level of pheromone sufficient to induce polarity. Next, the pheromone is rapidly removed (Figure 7A,B, top panels). For the NFB model, the system that started in the polarized state has a higher amount of $[eqn]$ at the time the pheromone concentration is removed (Figure 7C, top panel), and polarity is lost (Figure 7A, middle panel). Depending on the length of time before pheromone is removed, the system that is started in the unpolarized state will maintain polarity (Figure 7A, bottom panel) or, possibly, lose polarity at a slower rate. However, for the IFF model, the amount of $[eqn]$ is the same in both scenarios (Figure 7D, top panel) and polarity is lost at the same rate (Figure 7B, middle and bottom panels). If we vary the duration of pheromone exposure, for the NFB model the uniform initial condition requires a longer signal duration than the polarized initial condition for polarity to be lost following removal of pheromone (Figure 7C, bottom panel), whereas the IFF polarity is always lost for the uniform initial condition (Figure 7D, bottom panel).

4. Discussion

Cdc42 is a member of the Rho family of GTPases. A primary role for these signaling molecules is the regulation of the actin cytoskeleton. Mathematical models have been used to investigate many aspects of Rho GTPase signaling. The signaling circuit shown in Figure 2A not only underlies polarity establishment [16,40,41,42,43,44,45] but also produces waves and other patterns of signaling activity observed in cellular functions including cell division, cell migration, and phagocytosis [46,47,48,49,50]. The different behaviors arise by varying parameters in the core circuit that governs the strength and timing of feedback regulation. While these systems are typically modeled as mass-conserved reaction diffusion equations, intrinsic fluctuations arising from the stochastic nature of biochemical reactions can significantly impact their behavior. Therefore, understanding the mechanisms that regulate Rho GTPase signaling requires models that can capture multiple spatial and temporal scales.

4.1. Multiscale Models of the Polarity Circuit

Our modeling strategy was to start with a coarse-grained ODE model for the core polarity circuit that includes a known positive feedback loop in which Cdc42 recruits its own activator. To this core circuit, we added a negative feedback loop in which Cdc42 activity leads to the activation of a negative regulator. A key feature of the expanded model is that the feedback loops act on different time scales: rapid positive feedback and slow negative feedback. The computational and analytical tractability of the ODE model allowed us to demonstrate that it was sufficient to generate both bistable behavior and the rapid dismantling of the polarity site following the removal of the pheromone. While the ODE model provided evidence to support the role of negative feedback in dismantling the polarity site, it did not consider the spatial aspects of polarity establishment. Therefore, we turned to reaction–diffusion equations. Initially, to construct an RDE model, we simply added diffusion to our ODE model. In the absence of negative feedback this model reduces to the model studied by Chiou et al. [39] to investigate competition and coexistence between polarity sites. The model equations for the RDE model were simple enough that a qualitative understanding of the system’s behavior could be achieved through a combination of simulations and bifurcation analysis. Our analysis demonstrated that the qualitative properties of the system were maintained when spatial variations in concentrations were considered.

We then extended the model to represent the polarity circuit more faithfully by including more of the known chemical reactions in this system. Our extended model is essentially the one considered by Howell et al. [51] to study oscillations in the polarity site. While the extended model again qualitatively captured our experimental observations, the stochastic noise inherent in biochemical reaction and diffusion was not captured by the model. Previously, we used particle-based simulations of the polarity circuit to study the effects of intrinsic fluctuations on the performance of the system. Therefore, we converted the extend RDE model to a particle-based model using the method outlined in [52]. Again, the qualitative features of the system were preserved, and polarization occurred on time scales consistent with our experimental observations.

4.2. Negative Feedback Versus Incoherent Feedforward Regulation

Our results do not distinguish if pathway inhibition occurs through negative feedback or incoherent feedforward regulation. The biochemical network responsible for polarity establishment during mating also underlies polarity establishment during budding. There are known negative feedback loops that regulate polarization during budding and, therefore, are also likely to operate during the mating response. These include the reduction in activity of the GEF Cdc24 through feedback phosphorylation [53] and the recruitment of GAPs to the polarity site [54,55,56]. While there is no known incoherent feedforward loop that acts at the level of the polarity circuit, one does exist at the level of transcription regulation [36]. MAPK activity leads to activation of the transcriptional activator Ste12. At the same time, MAPK kinase activity leads to the activation of Far1, a multifunction protein. One function of Far1 is to increase the rate at which Ste12 is degraded [57], thereby forming an incoherent feedforward loop. Far1 is induced by Ste12 and also recruits the GEF Cdc24 to the polarity site. Therefore, the degradation of Ste12 could indirectly have a negative impact on polarity. While our experimental investigations were not sufficient to distinguish these potential mechanisms, recent computational and theoretical investigations have suggested methods for accomplishing this task in the case of ODE models [58]. Therefore, we investigated if a similar strategy could be employed with our models. Rahi et al. were interested in systems that adapt over time (i.e., show a transient response to a sustained signal). There are three key differences between the polarity circuit and the signaling networks considered by Rahi et al. First, the polarity circuit does not show an adaptive response; second, the polarity circuit is a spatially extended system; and third, the polarity circuit is bistable. We took advantage of the third property, bistability, to design an experiment that could, in theory, distinguish the two mechanisms. In the experiment, the system is prepared with two different starting conditions (polarized and unpolarized) and no pheromone present. In the case of NFB, the amount of negative regulator present in the system depends on the starting condition because the strength of negative regulation is tied to the concentration of active Cdc42. In contrast, for IFF regulation, the amount of the negative regulator is the same in both scenarios. Next, a pheromone concentration sufficient to establish polarity is added to both initial conditions. For NFB regulation, the negative regulator starts at an elevated level for the pre-polarized system, while it starts at a low level for the system in the unpolarized state. Therefore, if pheromone is removed soon after polarity is established, the system started from a pre-polarized state will immediately lose polarity. However, the system starting from the unpolarized state will either maintain polarity or lose it at a slower rate than the pre-polarized system. In contrast, for IFF regulation, the amount of the negative regulator will be the same for both starting conditions. Therefore, no difference will be seen when the pheromone is removed.

4.3. Potential Benefits of Negative Regulation of a Bistable Circuit

Bistability may provide important benefits during the mating response. Because yeast cells are not motile, they must undergo polarized growth to reach a distant mating partner. This chemotropic growth entails considerable energetic costs and presumably will be undertaken only when there is a high probability of successful mating. Monostable systems that contain strong positive feedback loops are characterized by sigmoidal or switch-like dose–response curves. Therefore, in the region where the polarity circuit transitions from an unpolarized to a polarized state, the system is highly sensitive to noise. Consequently, small fluctuations in the pheromone concentration can cause cells to gain or lose polarity. In contrast, bistable systems are more robust to perturbations and provide a mechanism for guarding against loss of polarity once polarity is established. Additionally, bistable regulation of polarity has also been suggested to improve sensitivity to shallow noisy gradients [59,60].

Our analysis also revealed the existence of negative regulation acting over time scales that are longer than the positive feedback loop driving polarity. Yeast cells often polarize in directions that are not aligned with the pheromone gradient and then reorient the polarity site as they grow. Because positive feedback tends to reinforce the polarity site, negative feedback is likely required to move the site. Indeed, previous studies have shown that vesicle delivery represents a delayed negative feedback loop by locally diluting Cdc42 concentrations [38,61,62] and is required for polarity patch movement. In another study that combined experimental and modeling investigations, it was demonstrated that negative feedback generates oscillatory behavior of the polarity site and increases the robustness of the polarity circuit to fluctuations in molecular abundances [51]. Our results suggest an additional physiological benefit of delayed negative regulation may be to rapidly dismantle the polarity site if the signal is abruptly lost when a competing suitor successfully mates with a nearby potential partner.

In summary, we applied a multiscale modeling approach to investigate feedback regulation in the yeast polarity circuit. Our strategy was to start with a deliberately simplified ODE model to gain insight into the system-level properties of the polarity circuit. Guided by our analysis of this model, we developed a series of spatiotemporal models of increasing biological detail and spatial resolution. While we focused on the yeast polarity circuit, we believe a similar approach can be applied to other intracellular signaling systems.