Models in biology: lessons from modeling regulation of the eukaryotic cell cycle

Jeremy Gunawardena [1] opened this series of review articles (on ‘models in biology’) with a quote from Nobel Laureate James Black, who described mathematical models as "accurate descriptions of our pathetic thinking about nature". Black, who made skilful use of mathematical models on the road to discovering antagonists of β-adrenergic receptors, was certainly not equating mathematical models with ‘pathetic thinking’ (as some people believe) but rather suggesting that mathematical models could be useful in turning our provisional ideas about molecular biology into real knowledge about living cells. Black saw three essential roles for models in biology: to expose assumptions, to define expectations, and to devise new tests. Our purpose in this review is to illustrate these principles by our experiences in modeling the control of mitosis in eukaryotic cells. We structure our review around subheadings derived from some cogent observations by Gunawardena (slightly paraphrased) about the art of building mathematical models in biology [1].

We focus in this review on the molecular mechanisms controlling entry into and exit from mitosis (M phase) in the natural cell cycles of frog embryos and fission yeast cells and in the artificial mitotic cycles exhibited by frog egg extracts. Because our focus is on certain principles of model building rather than a comprehensive description of cell cycle controls, we will not discuss the mechanisms controlling commitment of cells to DNA replication (S phase). Following an ancient tradition of cell biologists, we refer to G1, S and G2 phases of the cell cycle simply as ‘interphase’.

To describe the role of mathematical modeling in this endeavor, we start by summarizing the relevant knowledge — in 1990 — of the physiology and biochemistry of M-phase controls in frog eggs and fission yeast cells. Our description is based loosely on review articles of the time by Murray and Kirschner [2] and by Nurse [3].

Immature frog oocytes, which are arrested in G2 phase of meiosis I, can be induced by progesterone to proceed through meiosis I and arrest again in metaphase of meiosis II. These mature oocytes are ready for fertilization, which induces entry into a sequence of 12 rapid, synchronous, mitotic cycles culminating in the midblastula transition. Masui and Markert [4] demonstrated that the progesterone signal can be bypassed by injecting immature oocytes with cytoplasm from mature oocytes, which (they opined) contained a ‘maturation promoting factor’ (MPF) that triggered the activation of a pool of inactive MPF in the immature oocyte. Later experiments by Kirschner and colleagues [5] showed that ‘MPF activity’ is elevated whenever frog eggs are in M phase (meiosis in oocytes or mitosis in fertilized eggs), and the acronym was reinterpreted as ‘M-phase promoting factor’.

Using frog egg extracts, which recapitulate in vitro many features of MPF activation and inactivation, Murray and Kirschner, among others, showed biochemically that the rate-limiting step for MPF activation in extracts is the synthesis of a single protein, cyclin B, and that MPF inactivation, as the extract returns to interphase, is associated with cyclin B proteolysis by a ubiquitin-dependent pathway [6, 7]. (The anaphase promoting complex, APC, is the ubiquitin-conjugating enzyme that promotes cyclin degradation as cells exit mitosis.) In addition, MPF accumulates in cells during interphase in an inactive form (preMPF); then, as cells enter mitosis, preMPF is converted to active MPF in an autocatalytic (self-amplifying) process. These discoveries led to a picture of M-phase control in frog embryos similar to the diagram in Fig. 1a.

Meanwhile, Nurse and colleagues, studying cell growth and division in fission yeast by genetic tools, had discovered three genes involved in regulating entry into mitosis: cdc2, cdc25 and wee1. The genetic evidence suggested that the cdc2 gene product, the protein Cdc2, can be either ‘inactive’ (in interphase) or ‘active’ (in mitosis). Cdc25, the product of the cdc25 gene, appeared to promote the activation of Cdc2, and Wee1 to inhibit the activation of Cdc2. Significantly, Nurse found that the length of rod-shaped fission yeast cells when they enter mitosis is sensitively dependent on the dosage of the wee1 and cdc25 genes. Deletion of wee1 makes cells smaller, whereas overexpression of wee1 makes cells larger. Deletion and overexpression of cdc25 have the opposite effects. These facts are consistent with the notion that Wee1 is an inhibitor of entry into mitosis and Cdc25 is an activator. Upon cloning these genes, Nurse and coworkers found that Cdc2 is a serine/threonine protein kinase and Wee1 is a tyrosine/threonine protein kinase. The amino acid sequence of Cdc25 was unlike any other proteins known at that time. With characteristic caution, Nurse [3] diagrammed the genetic evidence as in Fig. 1b, although he speculated that Wee1 inhibited Cdc2 by phosphorylating it on neighboring threonine and tyrosine residues, and Cdc25 activated Cdc2 by promoting the removal of these phosphate groups.

The activity called MPF proved particularly difficult to purify biochemically, but this feat was eventually achieved by Lohka et al. [8], who showed that MPF activity was associated with co-purification of Cdc2 and cyclin B [9–11]. It was thought most likely that cyclin binds to Cdc2 to form a kinase-active heterodimer (a process called ‘stoichiometric’ activation), but the possibility that cyclin played a catalytic role in the activation of Cdc2 could not be ruled out.

Two other experiments are especially relevant to our story. Mark Solomon, working in Kirschner’s lab, was studying the activation of Cdc2 by cyclin B in frog egg extracts that lacked both cyclin synthesis and cyclin degradation [12]. Solomon supplemented these extracts with fixed concentrations of non-degradable cyclin B and measured the resultant activity of Cdc2 kinase (by its phosphorylation of a representative substrate, histone H1) as a function of cyclin B level. We shall call this a ‘signal-response’ curve: cyclin B level being the signal, and Cdc2-kinase activity being the response. Solomon found a distinct ‘cyclin threshold’ for Cdc2 activation (Fig. 1c). For cyclin concentrations less than the threshold, he observed only background levels of histone H1 phosphorylation, but as soon as the cyclin level was greater than the threshold, the extract abruptly exhibited a high rate of histone phosphorylation, which increased further with increasing cyclin concentration.

In a similar fashion, Marie-Anne Félix was studying Cdc2-induced degradation of cyclin in frog egg extracts in Karsenti’s lab [13]. To an interphase-arrested extract (no cyclin synthesis), she added a small amount of [³⁵S]-labeled cyclin B plus a measured amount of exogenously synthesized Cdc2:cyclin B heterodimers (i.e., MPF, measured in units of kinase activity per microliter of extract). She found that for MPF activity less than 6 U/μl, the [³⁵S]-labeled cyclin B was only slowly degraded, but for higher concentrations of MPF (≥12 U/μl), the radioactively labeled cyclin B was rapidly degraded, after a time delay of 15 minutes or longer (Fig. 1d).

In the next section, we describe three simple mathematical models of these interactions between Cdc2 and cyclin B. The strengths and weaknesses of these models, in light of the experimental results of the time, set the stage for later, more successful models.

In 1991 three competing models of the cyclin-MPF network were published by Norel and Agur [14], Goldbeter [15], and Tyson [16]. Basing their models loosely on the experimentalists’ informal diagrams sketched in Fig. 1a, b, all three modelers focused on certain features of the network that they thought to be most important (Fig. 2a–c, left), while neglecting other biochemical details known (or suspected) at the time. For instance, Norel and Agur assumed that cyclin drives the production of MPF in a catalytic manner and that MPF activates its own production (‘autocatalysis’). Norel and Agur also assumed that MPF activates cyclin degradation (by activating the APC, the enzyme ‘E’ in Fig. 2a, left), but that cyclin degradation is distinct from MPF destruction. Goldbeter also posited a catalytic function for cyclin activation of MPF, without the positive feedback loop whereby MPF activates itself. Cyclin degradation in Goldbeter’s model is a two-step process, whereby MPF first activates the APC (‘E_a’ in Fig. 2b, left), and then the APC drives cyclin degradation. In Tyson’s model, on the other hand, cyclin is a stoichiometric activator of Cdc2; i.e., cyclin binds with phosphorylated Cdc2 to form preMPF (Fig. 2c, left). preMPF is then converted into active MPF by the phosphatase activity of Cdc25. Because Cdc25 is activated by MPF, the conversion of preMPF to active MPF is a self-amplifying process in Tyson’s model. Tyson neglected the role of MPF in activating the APC; instead he assumed, without much experimental evidence, that only a phosphorylated form of cyclin was rapidly degraded.

Based on the reaction networks in Fig. 2 (left panels), the modelers converted their biochemical assumptions about the interactions of cyclin and MPF into differential equations for the rates of change of cyclin level and MPF activity. These mathematical equations describe the expected properties of the hypothetical mitotic control system given the assumptions behind each model, and these properties can then be compared with experimental observations to assess the reliability of the model. For instance, all three models exhibit spontaneous oscillations of MPF activity similar to MPF oscillations observed in frog egg extracts and to the mitotic cycles observed in early frog embryos. MPF oscillations in each model depend, of course, on the values assigned to the rate constants peculiar to each model, as illustrated in the middle panels of Fig. 2, which shows how the amplitude of MPF oscillations depends on the rate of synthesis of cyclin B (the rate constant k _s). Notice that, in each model, the rate of cyclin synthesis must be greater than a certain minimum value for MPF oscillations to appear, as observed by Murray and Kirschner [6].

The models differ, however, in the dependence of MPF activity on a fixed concentration of cyclin (right panels of Fig. 2), which is Solomon’s signal-response curve (Fig. 1c). The Norel/Agur model is inconsistent with Solomon’s curve because the model has no ‘upper’ steady state of high MPF activity. In Goldbeter’s model, MPF activity increases abruptly over a small range of cyclin concentrations (a phenomenon called ‘ultrasensitivity’ [17]). Goldbeter’s ultrasensitive signal-response curve is consistent with the abrupt rise of MPF activity observed by Solomon. However, Goldbeter’s curve has a flat plateau for cyclin concentrations greater than the threshold, whereas Solomon’s measurements show a steady rise in MPF activity for increasing cyclin concentration above the threshold. In contrast to Goldbeter’s sigmoidal signal-response curve, Tyson’s model predicts an S-shaped curve, which is also consistent with Solomon’s experimental results, provided one interprets Solomon’s ‘threshold’ with the lower turning point of the S-shaped curve. Although both Goldbeter and Tyson referred to Solomon’s paper, neither of them commented on the significance of Solomon’s ‘cyclin threshold for Cdc2 activation’ with respect to their computed ‘signal-response’ curves.

The 1991 models were simple, with two or three differential equations (one for each time-varying species) and a handful of rate constants (one or two for each biochemical reaction). However, they were of little value in providing a unified, accurate picture of what was then known about mitotic regulation in eukaryotes and in making reliable predictions. To overcome these problems would require a more careful accounting of the molecular biology and cell physiology of frog eggs and yeast cells.

In 1993, we proposed a model of mitotic controls in frog eggs and extracts [18] that combined the best features of Tyson’s and Goldbeter’s models, namely, stoichiometric binding of Cdc2 and cyclin B, positive feedback loops through Cdc25 and Wee1, and delayed activation by MPF of the APC. Compared to the 1991 models, our 1993 model (Fig. 3) included more reactions among additional biochemical species (e.g., the regulatory enzymes Wee1 and Cdc25) and, hence, more differential equations to solve and more rate constants to estimate. In return, the added complexity allowed us to simulate in great detail many characteristic features of autonomous mitotic cycles in frog eggs [18] and of size-controlled division cycles in fission yeast [19, 20]. For example, the 1993 model accounts nicely for the classic experiments of Solomon et al. [12] and Félix et al. [13], as can be seen by comparing the simulations in Fig. 3b, c with the experimental results in Fig. 1c, d.

The successful features of the 1993 model gave us the confidence to make some intrepid predictions.

Theoreticians, especially physicists, like simple, elegant models with few parameters (four parameters supposedly suffice to ‘fit an elephant’), but living matter, it would appear, is not so simple. The 1991 models of MPF regulation were parsimonious but inaccurate and ineffectual. Effective models of mitotic control (models that can account for a wide range of important physiological features and that can make trustworthy predictions) require dozen(s) of differential equations (one for each important protein and post-translationally modified form) and an attendant proliferation of unknown rate constants. Nonetheless, we were able to show [19] that the model in Fig. 3a, under quite reasonable conditions, can be reduced to just two differential equations and a few control parameters. This reduced model provides the desired accurate and parsimonious picture of many characteristic features of mitotic control in frog eggs and fission yeast cells. However, to dig deeper into the molecular foundations of cell cycle control, including regulation of DNA synthesis as well as mitosis, requires the creation of more detailed dynamical models, as illustrated by our later studies of fission yeast [23] and budding yeast [24].

If we are compelled to build detailed models, then we must eventually face the problem of how to estimate dozens of unknown parameter values (rate constants, binding constants, enzyme activities, and so on). There is nothing new or unexpected about this problem; every rate constant is estimated ultimately by fitting a mechanistic model to experimental data. The more complicated the mechanism and the longer the list of unknown rate constants, the more data will be needed to estimate them. To estimate a large number of parameters, we will need reliable data that probe every aspect of the mechanism, but we do not need great quantities of data. A few data points on a time course of cyclin accumulation or disappearance may be enough to estimate the rate constant of cyclin synthesis or degradation with enough accuracy for our purposes. Indeed, we found, in fitting the 1993 model to frog egg-and-extract data, that we could estimate the important rate constants in the model with some confidence. One might consider these rate constant values as additional predictions of the model, since they were not, at the time, measured by independent biochemical investigations. Over the next few years, Kumagai, Dunphy and other biochemists, employing clever manipulations of frog egg extracts, provided direct measurements of many of the rate constants estimated by the model [25]. Although these biochemists were unaware of our model and certainly not intending to confirm our parameter values, their measured values are very much in line with our estimates [26].

Another surprising development is the recent recognition that one of the MPF-counteracting phosphatases, PP2A:B55δ, is regulated during entry into and exit from mitosis [34, 35]. This phosphatase (which we shall refer to as ‘B55’) is inhibited by the small proteins Arpp19 and endosulfine-α (which we will refer to collectively as Ensa). Ensa, in turn, is activated by phosphorylation by Greatwall kinase (Gwl), and Gwl is activated by phosphorylation by MPF (Cdc2:cyclin B) [36]. Experimental evidence suggested that B55 opposes the MPF-catalyzed phosphorylations of Wee1 and Cdc25 [34, 35] and of the APC [37], and we proposed (for theoretical reasons) that B55 is the Cdc2-counteracting phosphatase for Gwl as well. These known and predicted roles of B55 create a multitude of new positive and negative feedback loops in the network (Fig. 4a) with interesting dynamical consequences [38].

Of special significance is the double-negative feedback loop between Gwl and B55, which (we propose) creates a bistable switch controlled by MPF [38]. (Let’s call it the BEG switch, for B55-Ensa-Gwl [39].) As a cell enters mitosis, B55 activity is high, and high MPF activity is required to flip the BEG switch to the state of low B55 activity. Once the switch is flipped, moderate or even low activity of MPF is enough to keep B55 activity off and APC activity on (Fig. 4b). For this reason (we propose), as cells degrade cyclin B during anaphase and telophase, they do not prematurely inactivate the APC and abort cyclin B degradation. There is experimental evidence that APC activity functions as an ultrasensitive switch [40]; we propose that, on closer examination, APC will behave as a bistable switch (Fig. 4c) [38, 37].

If we correctly understand the role of regulated B55 in entry into and exit from mitosis, then one might reasonably ask how we could have achieved such a good description of mitotic controls in frog eggs with the 1993 model, which neglected any regulation of the MPF-counteracting phosphatase. The answer seems to lie in our choice of parameter values. In the 1993 model, we introduced a Cdc2-counteracting phosphatase (called ‘PPase’ in that paper, which we now know to be B55), and assigned it a moderate activity needed to set the limits of the bistable region (Fig. 3b) created by the interactions of MPF with Wee1 and Cdc25. Meanwhile, APC activation was controlled by a time-delay mechanism through the ‘intermediary enzyme’ IE in Fig. 3a. By adjusting the parameter values in the ‘IE-APC’ part of the network we could get APC to turn on after a time delay consistent with Félix et al. [13] (Fig. 3c) and to stay on long enough to degrade most of the cyclin B during the transition from metaphase to interphase. Compared to the 1993 model, the new model [38], with regulated B55, is parametrized differently and gives a better description (we believe) of the course of events at mitosis in frog embryos.

It is worth mentioning that the BEG switch is not essential for entry into mitosis in mammalian cells [39, 41]. In this sense, our 1993 model of frog egg extracts is actually more appropriate as a model of mammalian cells than of frog eggs, for which it was intended. Or we might say it is a model of frog eggs with constitutively phosphorylated Gwl.

The BEG mechanism was discovered because it is essential for mitotic entry in frog egg extracts [36, 42]. If Gwl and/or Ensa are depleted from extracts, then nuclei are blocked in late interphase (with replicated DNA). Why, we might ask, is Ensa essential in embryos but not in somatic cells? Recall that immature oocytes must arrest in G2 phase of meiosis I, awaiting a hormonal signal to undergo maturation. In this state, they may accumulate very high levels of preMPF, but they must not enter meiotic M phase prematurely. Having a high level of PP2A:B55δ could be a fail-safe mechanism for stabilizing the arrested oocyte. But in this case, B55 activity must be down-regulated in meiosis I in order for the egg to proceed into meiosis II, which might be the role of Gwl and Ensa.

If mathematical modeling is such a nifty tool, why has it taken so long for molecular biologists to incorporate it into their armamentarium? In our experience, mathematical modeling of the cell cycle was at first politely ignored, then actively opposed, and lately grudgingly accepted. Perhaps this is to be expected of any new approach that is unusual and difficult to employ. But most people are comfortable now with the realization that molecular cell biology is a complex, interdisciplinary field that requires expertise from many quarters: genetics, biochemistry, nanotechnology, image processing, and bioinformatics, for example. Modern biological research is team-based, and, where appropriate, these teams need expertise in mathematical modeling of molecular regulatory networks. This is a specialized skill, much like X-ray crystallography, molecular genetics, proteomics or super-resolution microscopy. Life science departments need to acknowledge the vital role of models in biology and hire computational biologists who can bring this tool to bear on the important and difficult problems of today and tomorrow.