Dynamic modeling of yeast meiotic initiation
© Ray et al.; licensee BioMed Central Ltd. 2013
Received: 31 December 2012
Accepted: 17 April 2013
Published: 1 May 2013
Skip to main content
© Ray et al.; licensee BioMed Central Ltd. 2013
Received: 31 December 2012
Accepted: 17 April 2013
Published: 1 May 2013
Meiosis is the sexual reproduction process common to eukaryotes. The diploid yeast Saccharomyces cerevisiae undergoes meiosis in sporulation medium to form four haploid spores. Initiation of the process is tightly controlled by intricate networks of positive and negative feedback loops. Intriguingly, expression of early meiotic proteins occurs within a narrow time window. Further, sporulation efficiency is strikingly different for yeast strains with distinct mutations or genetic backgrounds. To investigate signal transduction pathways that regulate transient protein expression and sporulation efficiency, we develop a mathematical model using ordinary differential equations. The model describes early meiotic events, particularly feedback mechanisms at the system level and phosphorylation of signaling molecules for regulating protein activities.
The mathematical model is capable of simulating the orderly and transient dynamics of meiotic proteins including Ime1, the master regulator of meiotic initiation, and Ime2, a kinase encoded by an early gene. The model is validated by quantitative sporulation phenotypes of single-gene knockouts. Thus, we can use the model to make novel predictions on the cooperation between proteins in the signaling pathway. Virtual perturbations on feedback loops suggest that both positive and negative feedback loops are required to terminate expression of early meiotic proteins. Bifurcation analyses on feedback loops indicate that multiple feedback loops are coordinated to modulate sporulation efficiency. In particular, positive auto-regulation of Ime2 produces a bistable system with a normal meiotic state and a more efficient meiotic state.
By systematically scanning through feedback loops in the mathematical model, we demonstrate that, in yeast, the decisions to terminate protein expression and to sporulate at different efficiencies stem from feedback signals toward the master regulator Ime1 and the early meiotic protein Ime2. We argue that the architecture of meiotic initiation pathway generates a robust mechanism that assures a rapid and complete transition into meiosis. This type of systems-level regulation is a commonly used mechanism controlling developmental programs in yeast and other organisms. Our mathematical model uncovers key regulations that can be manipulated to enhance sporulation efficiency, an important first step in the development of new strategies for producing gametes with high quality and quantity.
The diploid yeast Saccharomyces cerevisiae undergoes mitosis in glucose medium. Upon transfer to acetate sporulation medium, cells commit to meiosis, a division process that produces four spores . Meiotic initiation involves a sequential activation of signaling molecules. Importantly, expression of these molecules occurs transiently within a short time window [2–8], suggesting that protein turnover and modification are under tight regulation. These short-lived signals are important for efficient entry and successful completion of meiosis . Further, interactions among these signaling molecules can lead to different levels of sporulation efficiency, as seen from yeast strains with distinct mutations or genetic background . Understanding how the transient signals are generated and trigger sporulation at different efficiency represents an important first step in the development of new strategies for producing gametes with high quality and quantity.
Inactivation of PKA under meiotic conditions leads to enhanced activity of Rim11, a kinase that phosphorylates Ime1 and Ume6 [2, 4, 7]. Phosphorylation stabilizes the Ume6-Ime1 complex, which is recruited to the promoters of early meiotic genes such as IME2 to activate their expression [17, 18]. Ime2 is a kinase and functions as a positive regulator for premeiotic DNA replication and nuclear division . Ime2 plays an important role in terminating expression of early meiotic genes through promoting proteolytic degradation of Ime1 [3, 9]. This negative regulation ensures that Ime1 is only expressed within a narrow time window. In addition, Ime2 transcriptionally activates its own expression via upstream activation sequences [17, 20]. This positive auto-regulation allows Ime2 to activate middle meiotic genes independent of Ime1. Further, mutual antagonism between Ime2 and G1 cyclins (Cln3/Cdc28) may be responsible for distinct response modes of meiotic genes to Ime1 levels .
NDT80 encodes a transcription factor that activates expression of mid-meiosis genes . The phosphorylated Ime1-Ume6 complex is insufficient to activate NDT80 due to the presence of a repressor, Sum1, on its promoter. Ime2 can activate expression of NDT80 by eliminating Sum1-mediated repression [23, 24]. Ime2 further phosphorylates Ndt80, allowing Ndt80 to promote its own expression by competing with Sum1 for binding on the promoter . In turn, Ndt80 is a transcriptional activator of IME2. Ndt80 boosts Ime2 activity during the middle stage of sporulation [26, 27], and premature transcription of NDT80 induces transcription of IME2.
Because of complex feedback regulation on the meiotic initiation pathway, mathematical modeling becomes an important tool to understand dynamic behaviors of signaling molecules and how their interactions ensure different degrees of sporulation efficiency. Feedback controls, which link the output of a circuit back to its input, are a key mechanism to stabilize cell-fate decision. Both experimental data and computational modeling suggest important roles of feedback loops in regulating mitotic entry and exit, cell growth, cell cycle, and pheromone pathways [28–33]. Negative feedback loops can generate oscillations or monotonic dynamics, while positive or double-negative feedback loops can produce bistability, i.e., having two coexisting stable steady states [34–36]. In the case of stronger positive or double-negative feedback loops, bistability can further lead to irreversibility, where a cell is locked in the post-transition state even after the stimulus disappears . Likewise, feedback loops may be responsible for transient dynamics of early meiotic proteins and different sporulation outcomes.
Boolean network models—a discrete method—have been developed to simulate the dynamics of meiotic initiation pathways. One study focuses on predicting sporulation efficiency upon gene deletions, and the other explains transient transcription of IME1 and IME2 by introducing two hypothetical repressors to shut down gene expression [37, 38]. Here, we develop an ordinary differential equation (ODE) model, a continuous method, to faithfully describe the nonlinear temporal dynamics of meiotic initiation pathways, incorporating numeric values of protein expression, kinetic rates, and time. The model depicts a signaling cascade of early meiotic proteins. Importantly, the model illustrates phosphorylation reactions and feedback loops that are crucial for directing the initiation of meiosis, based on the current knowledge available in the literature. The model correctly captures transient dynamics of meiotic proteins and accurately produces sporulation phenotypes of single-gene knockouts. We apply this model to investigate the contribution of feedback loops to transient behaviors of signaling molecules and sporulation efficiency.
Feedback loops in the model
Ime1 —| pSok2 —| Ime1
Ime1 — > pIme1 — > Ime2 —| Ime1
Ime2 — > Ime2
The constructed model is an abstract of real pathways, incorporating major players and events. The effects of other molecules are reflected indirectly in the model. For example, we assume PKA remains constant at a low level under meiotic conditions [39, 40]; the effects of upstream regulators of Ime1 (e.g., Cln3/Cdc28, Msn2/4, Snf1) are collectively represented by a general activation signal and a repression signal through pSok2; mutual inhibition between Ime2 and Cln3/Cdc28 and mutual activation between Ime2 and Ndt80 are both captured by positive auto-regulation of Ime2.
Parameters defined in the model
Synthesis, dimension =/hour
Synthesis rate of Ime1
Synthesis rate of Ime2
s ’ ime2
The maximum rate of auto-regulation-dependent Ime2 synthesis
Degradation, dimension =/hour
Degradation rate of Ime1
d ’ ime1
The maximum rate of Ime2-activated Ime1 degradation
Degradation rate of pIme1
Degradation rate of Ime2
Phosphorylation, dimension =/hour
Phosphorylation rate of Rim11
Dephosphorylation rate of Rim11
Phosphorylation rate of Ume6
Dephosphorylation rate of Ume6
Phosphorylation rate of Sok2
Dephosphorylation rate of Sok2
Phosphorylation rate of Ime1
Constant measuring half-maximum inhibition of Sok2 phosphorylation by Ime1
Constant measuring half-maximum inhibition of Ime1 synthesis by pSok2
Constant measuring half-maximum activation of Ime1 degradation by Ime2
Constant measuring half-maximum activation of Ime2 synthesis through auto-regulation
Constant measuring half-maximum activation of Ime2 degradation
Steady state value of variables using a mitotic initial condition and baseline parameter values
Steady state value
High-throughput screens of ~4,000 yeast deletion strains have identified 267 genes required for sporulation (sporulation-deficient genes) and 102 genes whose disruption enhances sporulation efficiency (sporulation-proficient genes) . Our mathematical model describes temporal dynamics of meiotic proteins encoded by five genes, among which RIM11, UME6, IME1, and IME2 are sporulation-deficient genes and SOK2 is a sporulation-proficient gene. Because the sporulation data are distinct from those used for model building, they can be applied for evaluating model performance. We virtually delete each of the five genes from the wild type model and simulate temporal dynamics of proteins in these knockout models (Additional file 2: Tables S1, S2, S3, S4 and S5, Figures S1, S2, S3, S4 and S5).
The overall pattern indicates that early meiotic proteins are sensitive to parameters that directly regulate their homeostasis. Levels of Rim11, pUme6, and pSok2 are mainly affected by phosphorylation and dephosphorylation. Rim11 and pSok2 are more sensitive to dephosphorylation than phosphorylation (u rim11 , u sok2 ), but the opposite is true for pUme6 (p ume6 ). The findings are consistent with the active forms of these proteins in meiosis. Synthesis and phosphorylation are most important to alter the Ime1 level (s ime1 , p ime1 ); both processes directly determine the gain or loss of Ime1. The level of pIme1 is primarily modulated by the synthesis of Ime1 and degradation of Ime2 (s ime1 , d ime2 ). Ime1 synthesis indirectly controls pIme1 through regulating Ime1; Ime2 degradation indirectly influences pIme1 through Ime2-activated Ime1 degradation. Parameters that control Ime2 auto-regulation and degradation have the greatest influence on Ime2 variations (s ’ ime2 , d ime2 ).
Feedback regulations are important for coordinated and transient behaviors of developmental systems [33, 37]. Ime1 and Ime2 exhibit orderly and transient expression during meiosis (Figure 2). These short-lived signals are critical for the successful completion of sporulation . We investigate how feedback loops affect dynamics of early meiotic proteins. A total of three feedback loops are described in the model: double-negative feedback between pSok2 and Ime1, negative feedback from Ime2 to Ime1, and auto-positive feedback of Ime2 (Table 1). We up-regulate, down-regulate, or delete each feedback loop through manipulating corresponding parameters. Protein dynamics are monitored in these in silico perturbation experiments.
Feedback regulations are known to control cell fate decision . In the context of yeast meiosis, feedback loops linking early proteins may be responsible for distinct sporulation efficiencies, traced by steady state levels of Ime2, the most downstream protein in the model. We perform bifurcation analyses on parameters governing feedback loops to determine which ones cause changes in the stability of Ime2 equilibrium.
PKA mediates phosphorylation of Sok2 and Rim11. To examine whether PKA is also a bifurcation parameter, we vary the phosphorylation rates of Sok2 and Rim11 simultaneously (Additional file 2: Figure S9). Two stable steady states, separated by an unstable steady state, are again observed. One stable state is the default equilibrium value of Ime2. When PKA activity is reduced, represented by lowering phosphorylation rates, the second stable state of higher Ime2 appears, corresponding to elevated sporulation efficiency. This result suggests that sporulation efficiency can also be improved by suppressing PKA activity.
Parameter c 1 is the half-maximum constant of Ime2 inhibition on Ime1. Changing c 1 results in two stable steady states, separated by a very short segment of unstable steady state (Figure 8C). Baseline value (c 1 = 0.01) leads to the default Ime2 equilibrium (Ime2 = 0.27). The inhibition from Ime2 to Ime1 decreases with the increase of c 1 , producing enhanced Ime2 level. This indicates that sporulation efficiency can be improved by repressing negative feedback from Ime2 to Ime1.
The Hill coefficient determines the switch-like behavior of Ime2 equilibrium. We find that the range of c 2 in which the system exhibits bistability is sensitive to the Hill coefficient (Figure 9). The system is monostable with coefficients of 1 or 3, since one value of c 2 corresponds to a single value of Ime2. Higher coefficients result in the transition from a monostable to a bistable system. A Hill coefficient of 7 expands the region of bistability across a broad range of parameter space, making the cell fate more robust with respect to perturbations in the feedback loop. This result indicates that cooperativity of Ime2 molecules is essential for producing bistable sporulation outcomes.
Precise regulation of a gene cascade in a coordinated manner is required for initiating a developmental program at the right time. This is often achieved through the activation of an upstream master regulator, which is controlled by multiple input signals and further regulates expression of downstream genes. Downstream genes, in turn, feed back to the regulator to modulate the entire pathway activity. The combinational nature of feedback loops ensures correct temporal dynamics of a developmental program [34–36].
The goal of this study is to understand and predict the effect of the control structure, i.e., feedback loops, on transient expression of early meiotic proteins and on distinct sporulation efficiencies observed in budding yeast. We construct a meiotic initiation pathway using an ODE-based model that includes regulation of Ime1, the master regulator, and five other early-meiotic proteins. We consider three feedback loops that control expression of these proteins: double-negative feedback between pSok2 and Ime1, negative feedback from Ime2 to Ime1, and auto-positive feedback of Ime2. In particular, Ime1 is controlled by an upstream inhibitor, Sok2, and a downstream inhibitor, Ime2.
The model is capable of simulating orderly and transient expression of meiotic proteins, without relying on putative repressors to shut down gene expression . The model is further validated by quantitative sporulation phenotypes of single-gene knockouts. We analyze the sensitivity of the model and find that proteins are sensitive to processes that directly regulate their levels. Subsequently, we perform in silico experiments on the model to understand the feedback mechanism on controlling transient protein expression and different sporulation efficiencies. The strength of mathematical models is that they serve as easily manipulatable systems for many perturbation experiments that are either extremely difficult or not tractable in a wet-lab setting.
The new insights gained from this study are two fold. First, we conclude that feedback loops play important roles in terminating expression of early meiotic proteins. Negative feedback from Ime2 to Ime1 is responsible for transient expression of both Ime1 and Ime2, in agreement with previous finding [3, 9, 37]. However, our study elucidates, for the first time, that the auto-positive feedback of Ime2 also ensures that Ime2 expression is confined to a narrow window. In our model, Ime2 responds in a graded mode to the Ime1 levels (Figures 5, 6 and 7), consistent with experimental observation that transcription of early meiotic genes is regulated by a gradient effect produced by Ime1 .
More importantly, the second new insight from exploration of the model is that feedback loops are responsible for tuning the efficiency of meiotic pathways. We perform bifurcation analyses on feedback loops using the equilibrium value of Ime2 as the pathway readout. We find that, by adjusting each of the two arms of mutual inhibition between pSok2 and Ime1, the system is able to move from a default meiotic state to a more efficient meiotic state. Similarly, by manipulating the strength of negative feedback loop from Ime2 to Ime1, the model readily produces a default meiotic state and a more efficient meiotic state. Auto-positive regulation of Ime2 is characterized by the Hill function with a high coefficient, providing a simple, reasonably accurate approximation for multiple regulations occurring on Ime2. This positive feedback generates a bistable pathway with two alternative stable steady states—the default meiotic state and a more efficient meiotic state. The robustness of bistability is sensitive to the Hill coefficient, indicating a strong cooperativity and nonlinearity in the response of Ime2 to the feedback. We propose that the combinational feedback regulation controls sporulation efficiency and guarantees that meiotic initiation proceeds in an accurate temporal scale.
Our mathematical model constitutes physical interactions of early meiotic proteins and provides mechanistic insights into ordered appearance of key regulators and sporulation efficiency. Such a model illustrates how different feedback regulations are integrated in the signaling pathway to cause changes in protein expression and meiotic outcome. The model is a reduced system of differential equations, including only Rim11, Ume6, Sok2, Ime1, and Ime2. Other proteins and/or links involved in meiotic initiation are traced indirectly. Validation using deletion mutants of meiotic genes suggests that major regulatory interactions have been captured. We demonstrate that the ordinary differential equation method can depict the most prominent features of signaling pathway during yeast meiotic initiation. Our mathematical model allows for uncovering key regulations that can lead to manipulation of the pathway to enhance sporulation efficiency. This represents an important first step in designing new strategies for producing gametes with high quality and quantity.
We develop a dynamic model to describe signaling pathways that operate during yeast meiotic initiation. Our study suggests that both positive and negative feedback loops control transient expression of early meiotic proteins, and multiple feedback loops regulate the efficiency of meiotic progression. Thus, yeast meiotic initiation is the consequence of systems-level feedback that leads cells into distinct sporulation states.
Phosphorylated or unphosphorylated Rim11, Ume6, and Sok2 are variables in Equations 1, 2 and 3. Because these three proteins exhibit uniform expression levels over the entire course of sporulation [2, 39, 40, 42], the total amount of phosphorylated and unphosphorylated forms of each protein is assumed to be constant (one unit). Phosphorylation and dephosphorylation are the only events described. PKA catalyzes phosphorylation of Rim11 and Sok2, although it is not included explicitly as a variable in the model. Phosphorylation of Sok2 is inhibited by Ime1, which is described using an inhibitory Hill function . Phosphorylation of Ume6 is mediated by Rim11 and modeled by mass action [4, 7].
Unphosphorylated and phosphorylated Ime1 are variables in Equations 4 and 5. Ime1 synthesis is inhibited by Sok2, as described using an inhibitory Hill function . Ime1 has a basal degradation rate, plus an Ime2-induced degradation [3, 9]. We assume that the negative feedback-dependent rate of Ime1 degradation is proportional to a Hill function. Phosphorylation of Ime1 is mediated by Rim11 and defined by mass action .
Equation 6 describes the rate of change of Ime2, the most downstream protein in the model. The synthesis of Ime2 depends on phosphorylated forms of Ime1 and Ume6 [17, 18]. Ime2 synthesis is further enhanced through positive auto-regulation [17, 20]. We assume that the auto-activation obeys cooperative kinetics, modeled by a Hill function with the coefficient of 5. Ime2 degradation occurs in a density-dependent manner.
The ODE model in the Systems Biology Markup Language format is provided as Additional file 3. The model is also available at the BioModels database (http://www.ebi.ac.uk/biomodels/, MODEL1303060000).
We consider a mitotic initial state. All variables except for pSok2 are set to zero because a very low level of meiosis-specific proteins could be detected during vegetative growth [3, 7, 8, 39, 40, 42]. Sok2 functions as a positive regulator of mitosis but as a negative regulator of meiosis [6, 12]. Thus, pSok2 is given the maximum level of one as the initial condition.
Parameters are either rate coefficients with a dimension of per hour or dimensionless constants (Table 2). Parameter values are estimated from the literature when they are available. When no data exist for parameters, values are manually explored over several orders of magnitude. Baseline values of parameters are determined by comparing model-generated output with experimental results in the literature [2–8, 25] and by constraining variable values in the range of zero and one. Outlined below is how we estimate parameter values.
Synthesis rates of genome-wide proteins have been calculated from experiments when yeast cells grow mitotically in glucose medium . Fold increase in the IME1 promoter activity has been measured when cells are transferred from glucose to sporulation medium [15, 44, 45]. Synthesis rate of Ime1 is estimated to be 10/hour during meiosis based on the above two measurements. Protein synthesis of sporulating cells has been monitored through deep sequencing of ribosome-protected mRNAs . We obtain synthesis rate of Ime2 based on its protein production relative to Ime1 and synthesis rate of Ime1.
The half-life of Ime1 is 0.5 hour in the presence of Ime2 and 1 hour in the absence of Ime2 . Thus, degradation rates of Ime1 in the presence and absence of Ime2 are approximated as 2/hour (ln2/0.5 hour ≈ 2/hour) and 1/hour (ln2/1 hour ≈ 1/hour), respectively, assuming that Ime1 degrades exponentially with time. Accordingly, both Ime2-independent and Ime2-activated Ime1 degradation rates are set to 1/hour. We assume the phosphorylation status does not influence Ime1 degradation; therefore, the same value is used for the degradation rate of pIme1. The half-life of Ime2 is much shorter, approximately 5 minutes . Accordingly, the degradation rate of Ime2 is calculated to be 8/hour.
Although no kinetic data exist for phosphorylation, the dephosphorylation rate is estimated to be higher than the phosphorylation rate for Rim11 and Sok2 because the unphosphorylated forms of proteins are active during meiosis [2, 6]. On the other hand, phosphorylated Ume6 is dominant under meiotic conditions; its rate constant for phosphorylation is higher than that of dephosphorylation .
The ODE model is implemented in MATLAB. Numerical simulation of the model is performed with MATLAB using a non-stiff solver ode45. Numerical results are confirmed with other non-stiff or stiff solvers: ode23, ode113, ode15s, ode23s, and ode23t.
We perform global sensitivity analyses using the software package SBML-SAT . A global sensitivity analysis explores the variation of model output to simultaneous perturbation of all parameters over a large range. The method of multi-parametric sensitivity analysis is utilized, which implements Latin Hypercube Sampling to randomly generate parameter values from a given range using uniform distributions. Parameter range is within one order of magnitude larger and smaller than baseline values; a total of 5,000 samplings are performed for each analysis. For each randomly generated parameter set, an objective function is computed by the sum of square errors between model outputs from random and baseline parameter sets. Each parameter set is classified as “acceptable” or “unacceptable” if the objective function value is smaller or larger than the average of all objective function values, respectively. The cumulative frequency is calculated for both acceptable and unacceptable cases. Parameter sensitivity is defined by the maximum distance of the two cumulative frequencies according to the Kolmogorov-Smirnov statistics. Therefore, sensitivity is in the range of 0 and 1; more important parameters have larger sensitivity values.
We use the software XPPAUT (http://www.math.pitt.edu/~bard/xpp/xpp.html) for bifurcation analysis. Baseline parameter values and initial conditions are applied. Numerically stable and unstable steady states of the ODE model are determined as a function of a bifurcation parameter.
The authors would like to acknowledge Leanne Whitmore for critical reading of the manuscript. This study was supported by the March of Dimes Basil O’Connor Starter Scholar Research Award #5-FY10-485 to PY.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.