A predictive model of gene expression reveals the role of regulatory motifs in the mating response of yeast

Cells use signaling pathways to receive and process information about their environment. These systems are nonlinear, relying on feedback and feedforward regulation to respond appropriately to changing environmental conditions. Mathematical models developed to describe signaling pathways often fail to show predictive power, because the models are not trained on data that probe the diverse time scales on which feedforward and feedback regulation operate. We addressed this limitation using microfluidics to expose cells to a broad range of dynamic environmental conditions. In particular, we focus on the well-characterized mating response pathway of S. cerevisiae (yeast). This pathway is activated by mating pheromone and initiates the transcriptional changes required for mating. Although much is known about the molecular components of the mating response pathway, less is known about how these components function as a dynamical system. Our experimental data revealed that pheromone-induced transcription persists following removal of pheromone and that long-term adaptation of the transcriptional response occurs when pheromone exposure is sustained. We developed a model of the regulatory network that captured both persistence and long-term adaptation of the mating response. We fit this model to experimental data using an evolutionary algorithm and used the parameterized model to predict scenarios for which it was not trained, including different temporal stimulus profiles and genetic perturbations to pathway components. Our model allowed us to establish the role of four regulatory motifs in coordinating pathway response to persistent and dynamic stimulation.


INTRODUCTION
Proper cellular function requires cells to respond appropriately to stimuli in their environment. Environmental cues, such as hormones and growth factors, are typically sensed by receptors on the cell surface and transmitted by intracellular signaling pathways. A key function of these pathways is to initiate the appropriate transcriptional program to respond to the environmental challenge. Mathematical modeling has helped to elucidate many of the design principles that regulate the spatiotemporal activity of signaling pathways and allow them to function reliably in changing environmental conditions (1). The ultimate test for these models is to predict pathway dynamics under conditions of time-dependent stimulation regimens and in the presence of genetic or pharmacological perturbations that disrupt the system in well-defined ways. While many models have reproduced qualitative features of signaling systems, their quantitative predictive power is often lacking. One reason for the lack of predictive power is that many previous studies have assessed cellular responses only to constant stimuli. However, signaling networks are nonlinear systems which typically have both positive and negative feedforward and feedback loops that operate on different time scales. Therefore, full characterization of these systems requires using time-dependent stimulus profiles that probe multiple time scales (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12).
We have performed such an analysis using the mating response of Saccharomyces cerevisiae (yeast). This response is activated when a mating-type specific pheromone binds to and activates a G-protein coupled receptor on a cell of opposite mating type. The signal is then propagated by a mitogen activated protein kinase validate the model and demonstrate its predictive power. First, we used the model to predict the behavior of mutations that selectively disrupt various signaling motifs in the pathway. Then we used the model to predict the transcriptional response of the system at a lower pheromone concentration. The result of our investigations is a fully validated model of transcriptional regulation that allows a quantitative characterization of the signaling motifs that regulate gene expression. We anticipate that our approach provides a template for a research strategy to characterize regulatory motifs inherent to many signaling pathways.

Adaptation and persistence in the mating response pathway
To determine the dynamics of the yeast mating response, we developed experimental tools that allow cells to be exposed to well-defined input signals of any specified temporal profile and a readout that faithfully tracks the dynamic response of the pathway. For controlling stimulus profiles, we employed a microfluidics system that is an adaptation of the "dial-a-wave" system developed by J. Hasty and colleagues (27). For tracking timedependent changes in pheromone-induced transcription in living cells, we placed a shortlived fluorescent reporter under the control of the pheromone responsive FUS1 promoter.
The fluorescent protein we used is fast maturing (~15 min) and through use of an Ndegron tag (YΔk) was engineered to have a half-life similar to its mRNA (~ 7 min) (28).
The short-lived reporter is essential in studies of temporal response dynamics, since it reveals transient response characteristics that are otherwise masked by accumulation of a long-lived reporter protein.
Initially, we exposed cells containing our short-lived fluorescence reporter to a constant stimulus of 50 nM pheromone for 10 hrs and monitored reporter fluorescence by imaging of cells in the microfluidic chamber ( Fig. 2A). For all the experimental results presented in this manuscript, we used cells lacking the protease Bar1 to remove the effect of pheromone degradation (29,30). We refer to this strain as wildtype hereafter. Under these conditions, the transcriptional response of wildtype cells reaches a maximum amplitude at 220 min, and then decreases for the remainder of the experiment (Fig. 2B).
In our next studies, we exposed cells with the short-lived reporter to pheromone pulses of different duration and again monitored reporter fluorescence (Fig. 2C).
Interestingly, reporter gene expression was significantly sustained following removal of pheromone for pulses of 90 min or less. We refer to this property as persistence and quantify it as the time from removal of pheromone to the time that signal drops below 2.5% of the maximum of each response cuve. The extent of persistence is negatively correlated with the duration of the stimulus pulse; as pulse length increases the persistence of the transcriptional response decreases (Fig. 2D). Another important observation is that the rate at which the fluorescent reporter decreases in time is independent of pulse duration (Fig. 2E) and the half-life associated with this rate (98 ± 9 min) is considerably longer than the half-lives (~7 min) of the reporter mRNA and protein ( Fig. 2E, green line) (28). Thus, new synthesis of transcripts and protein continues during the attenuation phase.
A simple explanation for the observed pathway persistence is that it represents a delay between receptor signaling and translation and maturation of the induced GFP reporter. To test this possibility, we developed a linear mathematical model of the response pathway that takes into account this delay (Supplementary Materials). Our analysis of the model revealed that a simple delay cannot account for the persistence in the transcriptional response (Fig. S1). In total, our preliminary investigations reveal that the pathway contains some form of "memory" that sustains new mRNA synthesis following removal of pheromone.
We next sought to determine at what level in the pathway the mechanisms for long term adaptation and persistent signaling occur. To determine if "long-term adaptation" relies on upstream pathway regulators of short-term desensitization, such as Sst2 or receptor endocytosis (31)(32)(33)(34), we investigated the dynamics of MAPK activity. We monitored MAP kinase dual phosphorylation, which is an indicator of activity, by Western blotting protein extracts of aliquots prepared from cells in the presence of 50 nM pheromone for a 10 hr time course. Fus3 activity remained constant after a transient increase and that of Kss1 increased throughout most of the time course and only slightly diminished toward the end of the experiment (Fig. 3A). These results demonstrate that the mechanism of long-term adaptation of transcriptional response does not involve upstream signaling events, but likely occurs at the level of transcriptional regulation.
We similarly monitored Fus3 and Kss1 kinase activity for a 90 min pulse of 50 nM pheromone by Western blot analysis for dual phosphorylation of the MAPKs. In this case aliquots of the culture were removed at indicated intervals during pheromone exposure and after removal of pheromone. Unlike gene expression, activity of the two MAPKs diminished rapidly once pheromone was removed (Fig. 3B), demonstrating that the mechanism for persistence also lies downstream of the MAPK signaling molecules.

Model for transcriptional regulation
To determine which elements of the pathway are critical for regulating the magnitude of the response, long-term adaptation and persistence, we developed a mechanistic model. We chose to include four established signaling motifs. We included an incoherent feedforward loop resulting from Far1-dependent Ste12 degradation as one potential mechanism for long-term adaptation (22) (Fig. 4A, motif 1). Next, we included positive feedback loops resulting from Ste12 auto-regulation and Ste12-dependent transcription of the MAPKs (Fig. 4A, motif 2), which we hypothesized could contribute to both the amplitude and persistence of the signaling response. Previously, we found that if Ste12 in the Dig/Ste12 complex degraded more slowly than free Ste12, rebinding of the Digs to Ste12 could act as a mechanism for adaptation ( Fig. 4A, motif 3) (35). We hypothesized that this motif also could contribute to persistent signaling, if the rate constant associated with rebinding was small. Finally, we included a negative feedback loop resulting from Ste12 induced synthesis of Far1 (Fig. 4A, motif 4), which we hypothesized might also contribute to long-term adaptation. These four motifs were included in the full model to capture persistent activation following stimulus removal and long-term adaptation of the transcriptional response (Fig. 4B). Importantly, the model also included synthesis and degradation of our transcriptional reporter. The abundance of this transcriptional reporter was the experimental output used to train the model. The model input was a piecewise-linear function that corresponded to the temporal pheromone stimulation profile. Full details of the mathematical model including the set of differential equations that describe the system are presented in Materials and Methods.
Below we describe the data sets used to train the model and present results for the model's performance. We used an evolutionary algorithm (Fig. 4C) to fit the model's 28 parameters. For each parameter, we determined a biologically relevant range from which the parameter values were selected (Supplementary Material). Each generation of the evolutionary algorithm had 500 individual parameter sets that underwent selection, crossover, and mutation. Over the course of 100 generations the total absolute error (TAE) between the experimental data and the simulations converged ( Fig. 4D and E).

Assessment of model performance
Signaling pathways represent nonlinear dynamical systems capable of responding on multiple different scales. Therefore, we reasoned to develop a predictive model for transcriptional regulation, it was critical to measure the system's response to time dependent pheromone concentrations with multiple different frequencies. To this end, in addition to the single pulse data described above (replotted in Figs. 5A-F), we collected data for periodic stimulation consisting of pulses of pheromone in which the on and off intervals were the same length. The on-off durations used were 45, 60, 75, 90, and 120 ( Fig. 5, G-K). We also included data for constant pheromone stimulation (replotted in Fig.   5L.) The model captured the varying durations of persistence after stimulus is removed ( Fig. 2D, gray curve) and long-term adaptation to stimulus under conditions of both periodic and constant stimulus (Fig. 5). The model was also capable of capturing the dynamics of the MAPK activation profiles in response to constant stimulation and to a 90 min pulse of 50 nM pheromone (Fig. S2).

Regulation of response to prolonged stimulus
To understand how response to constant stimulus is regulated, we perturbed motifs that are likely to affect the magnitude of the response and long-term adaptation. First, we used the model to investigate the role of Ste12 auto-regulation in determining the magnitude of the transcriptional response by eliminating motif 2. We found that some parameter sets found by the evolutionary algorithm predict a dampened response in the absence of Ste12 auto-regulation (green and brown curves) while other parameter sets predict a response similar to wildtype (purple curves) (Fig. 6A). To experimentally determine whether Ste12 auto-regulation in this system has a significant role in amplifying the response, we replaced the Ste12 promoter with that of the promoter of the scaffold protein Ste5 ( P STE5-STE12). We chose this promoter because it produces constitutive amounts of Ste12 similar to the basal amount from the endogenous promoter (Fig. S3) and is not subject to auto-regulation by Ste12 (36). In cells containing the P STE5-STE12 mutation, the overall transcriptional response was diminished, and long-term adaptation began ~50 min sooner than for wildtype cells (Fig. 6A, triangles). These findings indicate Ste12 auto-regulation is important for amplifying the response and affects the timing of adaptation. They also suggest that Ste12 autoregulation counterbalances the depletion of Ste12 promoted by Far1 (motif 1).
Next, we investigated the role of Far1 in long term adaptation. In the model, Far1 is involved in two signaling motifs. The first is the incoherent feedforward loop (motif 1) formed by MAPK activation of Far1, followed by active Far1 promoting the degradation of Ste12. The second (motif 4) is the negative feedback loop formed by Ste12-dependent expression of Far1. First, we used the model to predict the system's response when Far1 is eliminated, which blocks both Far1-dependent mechanisms of adaptation (motifs 1 and 4). We found that some parameter sets predict no long-term adaptation in the absence of Far1 (blue and brown curves); however, other parameter sets predict no difference from wildtype (purple curves) (Fig. 6B). We reasoned that the interaction between the Digs and Ste12 was responsible for long-term adaptation for those parameter sets still exhibiting adaptation in the absence of Far1 (35). In the model, we allowed for the possibility that complex formation with the Digs protects Ste12 from degradation. There are two consequences of this protective complex (35). First, it ensures a large pool of inactive Ste12 is maintained prior to pheromone stimulation. Second it provides for an adaptive response. The basis for adaptation is that following exposure to pheromone, the free Ste12 concentration transiently increases as Ste12 is released from the Digs, but eventually returns to its pre-stimulus level (35). To test whether eliminating protective binding has any effect on adaptation in our model, we set the degradation rate of Ste12 in the Dig/Ste12 complex equal to that of the degradation rate of free Ste12. With this change, long-term adaptation was lost when the model was run using parameter sets that predicted Far1 was not involved in adaptation (Fig. S4). To determine which of these mechanisms is responsible for regulating long-term adaptation, we examined the response of a far1∆ mutant strain. In this mutant the transcriptional output does not diminish over time (Fig. 6B, triangles) demonstrating that Far1-dependent degradation of Ste12 is the primary mechanism of long-term adaptation.
To further constrain model parameters, we retrained the model including experimental data for the P STE5-STE12 and far1Δ mutants. The resulting parameter sets better captured the responses of the pathway mutants than those used for predictions (compare  Fig. S6A and B). This demonstrates that including strategic pathway perturbations in the training data can improve ability to identify biologically relevant parameters.
Because elimination of Far1 disrupts both the incoherent feedforward and negative feedback motifs (motif 1 and 4, respectively), we used the model to test if negative feedback contributes to long-term adaptation. In the model, disruption of the incoherent feedforward loop is equivalent to eliminating Far1 since promoting degradation of Ste12 is the only effect of Far1 on transcriptional response. However, we can use the model to identify the role of negative feedback. When transcriptional induction of Far1 by Ste12 was eliminated (motif 4) some simulations predict a sustained response, but most simulations still show adaptation (Fig. 6E). These results suggest that the incoherent feedforward loop (motif 1) is the predominant mechanism for long-term adaptation of the transcriptional response. While the model did not require the induction of Far1 for long term transcriptional adaptation, it is likely this feedback is required for one of the other functions of Far1 in the mating response, such as gradient sensing or maintaining cell cycle arrest.

Regulation of response to dynamic stimulus
To examine motifs that could contribute to persistence and further test the model's predictive power, we measured the system's response to single 90-minute pulses of 50 nM pheromone in the presence of pathway mutants that perturb Ste12 autoregulation, binding to DNA, or binding to the Dig1 and Dig 2 repressors. First, we eliminated Ste12 autoregulation (motif 2) as before by using the P STE5-STE12 mutation. While there was a dampened response consistent with the response to constant stimulus, neither the simulations nor the mutant response had any appreciable effect on persistence (Fig. 7A).
Second, we examined a mutation to one of the three pheromone responsive elements (PREs) within the FUS1 promoter that drives transcription of the GFP reporter. Ste12 has been reported to bind at a synthetic promoter having the same PRE mutation with only 30% of the affinity that it binds to a synthetic promoter with the wildtype sequence (37).
This PRE mutation (PRE*-GFP) significantly reduced the maximal amplitude of the transcriptional response (Fig. 7B, triangles). Using the best 10% of parameter sets found from fitting to the wildtype, far1Δ, and P STE5-STE12 data, we predicted the response of the PRE mutant by increasing the apparent dissociation constant by a factor of 3.33. The resulting parameter sets accurately predict the response of the PRE*-GFP mutant (Fig.   7B).
Finally, we examined the effect deleting the Dig1 and Dig2 repressors, which causes Ste12 to be constitutively active. This deletion mutant showed high basal expression and a slight increase in expression following pheromone induction (Fig. 7C,   triangles). We predicted the response of the dig1Δdig2Δ mutant by setting the total Dig concentration to zero. The model predicted a wide range of responses in the absence of the Dig1 and Dig2 transcriptional repressors. Interestingly, the results could be clustered into three groups (Fig. 7C, colored curves). The model predictions that showed high basal transcriptional response (Fig. 7C, green curves) result from parameter sets in which the degradation of Ste12 in complex with the Digs (kdegS12D) is similar to that of the degradation rate of free Ste12 (kdegS12) (Fig. 7E). In this case the total amount of Ste12 is the same in the dig1Δdig2Δ mutant and wildtype reference. Removing the Dig repressors generates more active Ste12 prior to pheromone stimulation, and, therefore, higher levels of the reporter in the mutant. For model predictions in which the prestimulation level of the fluorescent reporter does not increase significantly compared to the wildtype reference ( Fig. 7C, cyan curves), removing the Digs had two effects. For these parameter sets, the degradation rate of Ste12 in complex with the Digs is reduced (Fig. 7E). That is, the Dig repressors provide protective binding. Removing the Digs exposes Ste12 for degradation, but also activates Ste12. When these two opposing effects are balanced, the pre-stimulation level of active Ste12 in the dig1Δdig2Δ mutant is similar to that of wildtype, and, therefore, the expression level of the reporter does not significantly increase. The parameter sets that fit the experimental data best had intermediate degradation rates for Ste12 in complex with the Digs (Fig. 7E and Fig. 7C, brown curves). These results are consistent with our previous analysis of Ste12 dynamics that demonstrated that the Dig repressors provide some degree of protective binding (35).
Another observation consistent with previous experimental observations is that parameter sets best fitting the experimental results for the dig1Δdig2Δ mutant (Fig. 7C, brown curves) predict that the rate at which active Far1 is degraded (kdegPF1) is less than that for inactive Far1 (kdegF1) (Fig. 7E) (38). Additionally, other parameters that affect Far1-dependent degradation of Ste12 including the rate of Far1 dephosphorylation Because in the model, the only mechanism for transcriptional induction is dissociation of Ste12 from the Dig repressors, the model is not able to capture the slight pheromone-dependent induction seen in the dig1Δdig2Δ strain. This induction may result from pheromone-induced degradation of the transcription factor Tec1, a known binding partner of Ste12 (39). The slight pheromone-dependent induction in the dig1Δdig2Δ strain exhibits prolonged maximal expression after a 90-minute pulse of stimulus (72 min persistence) compared to wildtype (43 min persistence). To further investigate how the transcriptional repressors contribute to persistence, we perturbed motif 3 by increasing the rebinding rate of Ste12 to the Digs in the model by 5-fold. In doing so, the average persistence of the simulations decreased from 25 min to 13 min (Fig. 7D). This result combined with the prolonged persistence when the transcriptional repressors are deleted suggests that slow rebinding of the transcriptional repressors are a primary factor in the persistent transcriptional response following stimulus removal.

Prediction of different stimulation profiles
To further test the model's predictive power, we measured the response of cells exposed to periodic stimulation at the same frequencies as shown in Fig. 5, but at a lower pheromone concentration (Fig. 8, triangles). In response to 10 nM of constant pheromone, the fluorescent reporter achieves the same maximum amplitude as the 50 nM case but takes 25 min longer to reach its half maximum amplitude (Fig. 8A). For short pulses of stimulus (Fig. 8B), the amplitude of the response to 10 nM is considerably lower than that to 50 nM for all pulses. However, for longer pulses (Fig. 8C) there is less of a difference in the amplitude between the two doses. To simulate the lower pheromone dose, the only modification we made to the model was to adjust the slope of the input signal to match the slower production rate of the fluorescent reporter measured at 10 nM constant pheromone. Using this adjustment to the input stimulus, all of the parameter sets that fit the 50 nM data accurately predicted the response to sustained and pulsed 10 nM pheromone (Fig. 8, blue curves). This performance demonstrates that this model is capable of capturing behaviors at different doses of stimulus despite only being trained on a single dose.

DISCUSSION
A common way for cells to respond to changes in their environment is by regulating gene expression. Because environmental conditions are dynamic and can show significant variability, gene expression needs to be tightly regulated by the signaling pathways used by cells to monitor their surroundings. This regulation, which typically takes the form of feedback and feedforward loops, makes gene regulation an inherently non-linear process.
Therefore, predicting the response of these systems is not possible without the aid of mathematical models. Developing predictive models is challenging for two reasons: 1) these models tend to contain many parameters that are not directly measurable and therefore must be estimated from experimental data and 2) these systems operate on multiple time scales and, therefore, experimental data sets used to train the models must capture the relevant time scales. To overcome both these obstacles, we developed a research strategy that involved exposing yeast cells to single and periodic pulses of mating pheromone. By varying both the duration and frequency of the pulses, we ensured that the regulatory network that controls gene expression during the mating response was probed on the relevant time scale and with sufficient temporal resolution to accurately perform parameter estimation. This systematic analysis revealed novel features of the pathway and allow us to develop a mathematical model with predictive power.
Our analysis led us to the discovery of memory in the yeast mating response.
Specifically, we discovered that transcriptional regulation was sustained following removal of pheromone, a property of the system that we refer to as "persistence". Our model revealed that this persistence was not due to positive autoregulation of Ste12 but rather involves slow rebinding of the transcriptional repressors to Ste12. Persistent signaling may represent an important design feature of the pheromone response pathway. Yeast mating takes place in noisy environments where pheromone levels are expected to fluctuate. Preparing for mating takes a significant fraction of the cell's resources. Therefore, once the decision has been made to commit to the mating, it is important that the cell not "give up" if there is a transient loss of the pheromone signal.
Persistent signaling provides a mechanism to guard against this situation. Conversely, it is also important that a cell not remain committed to mating indefinitely. This might explain why persistent gene expression does not rely on positive feedback, which is capable of generating an irreversible switch.
In the presence of sustained pheromone signals, it is probably beneficial for yeast cells not to remain growth arrested when mating is unlikely to be successful. Such adaptative behavior in the mating response has been observed previously (22,40).
However, these studies were done in the presence of the protease Bar1, which degrades pheromone, thus making it difficult to identify the predominate mechanism that underlies transient signaling. Interestingly, our results revealed that in the absence of Bar1, pheromone-induced gene expression is transient, whereas MAPK signaling is sustained.
Our model predicted two mechanisms could underlie this long-term adaptation, an incoherent feedforward loop or protective binding. Further experiments revealed that the incoherent feedforward loop involving Far1-dependent degradation of Ste12 accounts for most of the long-term adaptation.
Our approach that combines mathematical modeling with experiments designed to probe cellular response pathways over multiple time scales provides a general framework for investigating gene regulatory motifs. First, we used experiments to narrow the portion of the pathway responsible for the dynamic properties of interest. In our case, these preliminary investigations revealed that both persistent signaling and long-term adaptation occurred at the level of gene regulation and did not involve upstream signaling components. Next, we developed a model incorporating known regulatory mechanisms and narrowed parameter ranges to physiologically relevant values. We then performed parameter estimation using an evolutionary algorithm applied to training data sets spanning multiple timescales. The use of time-dependent stimuli covering multiple time scales was essential for building a predictive model. When a subset of the data was used, model parameters were significantly less constrained, and the model's predictive power was reduced. Additionally, training on data sets spanning multiple timescales revealed the differences in timing of the signaling motifs. For example, the rebinding of the transcriptional repressors and the incoherent feedforward operate on different timescales, leading to decreased persistence after longer pulses of stimulus.
We also note that successful model building is an iterative process. For example, when fit only using wildtype data the model found two mechanisms of long-term adaptation were consistent with the data. The model also predicted positive feedback contributed to amplifying the signal but showed significant variability in the predicted strength of this feedback. Strategic experiments using targeted mutants were then able to identify the true mechanism of long-term transcriptional adaptation and quantify the role of positive feedback. Including these results in the training data sets, further constrained parameter values and allowed the model to accurately predict the system's behavior for lower pheromone concentrations and additional genetic perturbations.
Because gene editing and quantitative experimental approaches are becoming increasingly more feasible in other cell types, including mammalian cells, we believe our approach can be adapted to these systems. For example, such studies could reveal important information about the dynamics of MAPK signaling pathways dysregulated in diseases, including cancer, and ultimately suggest treatments for restoring proper function. Table 1 lists plasmids used in this study. Those that have been described previously are listed with the corresponding reference. Standard recombinant DNA procedures were used for construction of those plasmids described below (41). Table 2 lists the sequence of oligonucleotides used for PCR fragment amplification, mutagenesis, and DNA sequence confirmation involved in plasmid and strain constructions.

Cell Extract Preparation and immunoblotting
The following procedure was performed to determine the phosphorylation state The fluorescent signal is quantified in the chamber as the height of the input channels is alternated. Typically, it takes 2-5 min for the dye to equilibrate inside the chamber after switching the channels. This is much faster than the timescale of transcriptional response in the mating pathway. During the experiments, switching is automated using a step motor. Detailed microfluidics can be found in a methods review (51).

Microscopy
All experiments were performed in a microfluidic device using cell culturing methods as

Image analysis
To improve the quality of segmentation images were edited in ImageJ by first subtracting the background, using the unsharp mask filer with a radius of 2.0 and a mask weight of 0.5, using the Gaussian blur filter with a radius of 3.0, and finally subtracting the background again. These images were then used to perform image segmentation using and G2/M phases of the cell cycle when cells were exposed to 50 nM pheromone. This comparison revealed that cells in all phases of the cell cycle respond to pheromone by inducing transcription but those in S-phase respond more slowly than those in G1 and G2/M (Fig. S7).

Model Development
Our experimental investigations revealed that the mechanisms responsible for persistent gene expression following removal of pheromone must occur downstream of MAP kinase activity. Therefore, we do not explicitly consider (1.1) where !"#$% is the constitutive synthesis rate of Fus3, 123$% is the MAPK degradation degradation. The equations that govern the concentration of free Ste12, [ 12], and Ste12-Digs heterodimer, : 12 453! ; are given by: represents the total Dig1/2 concentration, which is assumed to remain constant, 9( is the association rate constant, krsd2 is the dissociation rate constant in the absence of pheromone and 90 [ ] is MAPK dependent increase in dissociation rate.
Since degradation of Ste12 is dependent on active Far1, our model includes Far1 dynamics. The equations that govern the concentrations of active Far1, [ 1], and inactive Far1, [ 1] are given by: (1.5) where !"#$( is the constitutive synthesis rate of Far1, /% is the dephosphorylation rate where !"#>$-is the constitutive synthesis rate of GFP, 123>$-is the degradation rate of GFP, which is known as we used a short lived GFP with a well-established half-life, and the term 9= [ 12] )* >$-+ [ 12] )* ⁄ models synthesis of GFP due to Ste12 dependent gene transcription.

Modeling pheromone signal
In our model the pheromone signal comes in at the level of the MAPK, Fus3, which is activated in a pheromone dependent manner. We model upstream activation of the pathway as a piecewise linear function. We specify a slope ?# that describes the rate at which signal activity increase following pheromone exposure and assume the signaling turns of instantaneously following removal of pheromone. The maximum input signal activity is kp1. For pulse trains of period p (on+off phase), the maximum input signal ( 89@ ) achieved during the simulation ( 89@ ) is: Then we use the following piecewise linear equation to describe the temporal signal profile for periodic stimulus: where is the time and is the pulse number defined by floor( / ).
For single pulses of stimulus, the signal is described by the following equation:

Parameter estimation
To perform parameter estimation, we used an evolutionary algorithm. Evolutionary algorithms have the goal of optimizing a solution to a problem and are inspired by elements of biological evolution including recombination, mutation, and selection. In our application we aim to optimize the fit of our model to experimental data by finding parameter sets that minimize the error between the experimental and simulated data. We implemented the algorithm using DEAP (Distributed Evolutionary Algorithms in Python) which is a user friendly framework for building and executing evolutionary algorithms (53).
The evolutionary algorithm is broken down into three main functions: simulation, scoring, and evolution. The algorithm is initiated by selecting parameter sets from specified uniform random distributions. In the simulation function, these parameters are used to simulate the model. In the scoring function, the difference between the simulation and the experimental data is the quantified using the mean absolute error. In the evolution function, the best parameters are chosen through a tournament. Then those parameters go through mating and mutation. Mating was simulated using a two-point crossover function and mutation was simulated using a polynomial bounded mutation function. The resulting parameter sets are then returned to the simulation function and the process starts over again. The algorithm continues to optimize parameter sets for a specified number of generations.
We performed all model fitting with 100 generations of 500 individuals because these numbers were typically sufficient for convergence of the score function. The mating and mutation functions were chosen because they worked best with synthetic data sets, specifically to optimize models with hill functions. Similarly, the hyperparameters (mutation rate, crossover rate, and tournament size) were selected based on their efficiency in fitting the synthetic data set best.

Scaling experimental data
It was necessary to scale the experimental data sets to account for differences in experimental conditions, such as light source and intensity. We first normalized the wildtype time series for constant 50 nM ⍺-factor to have a maximum of 1. Then all data sets generated using 50 nM ⍺-factor were scaled to align with the wildtype response during the initial on phase of pheromone exposure. For example, the choice of scaling factor for a 45 min pulse experiment was based aligning the first 45 minutes of this times series with the response for the constant pheromone case.
For the pathway mutants, a series of constant stimulation experiments was done in quick succession using the same microscope settings. In this experiment data were collected for wildtype, far1Δ, dig1Δdig2Δ, PRE*-GFP, and P STE5-STE12 strains. The wildtype data was multiplied by a scaling factor to align with the data used to train the model. The mutant responses were then multiplied by this scaling factor to correctly adjust their starting and maximal values. The response of the each of the mutants to a 90-miunte pulse of stimulus was then multiplied by a strain specific scaling factor to match the corresponding scaled constant response.

Prediction response to low pheromone dose
To predict the low dose data set, we scaled the signal input signal on rate ( ?# ) by 0.3 based on an estimation of the difference in the slopes for the temporal response to 50 nM and 10 nM constant pheromone. Using the parameter sets corresponding to the best fits to the 50 nM data but changing ?# , we successfully simulated the 10 nM data ( Fig. 8, blue curves).

Predicting response of pathway perturbations
The model was used to predict the response of four mutations to the pathway,