Mathematical Model Explaining the Role of CDC6 in the Diauxic Growth of CDK1 Activity during the M-Phase of the Cell Cycle.

In this paper we propose a role for the CDC 6 protein in the entry of cells into mitosis. This has not been considered in the literature so far. Recent experiments suggest that CDC 6 , upon entry into mitosis, inhibits the appearance of active CDK 1 and cyclin B complexes. This paper proposes a mathematical model which incorporates the dynamics of kinase CDK 1 , its regulatory protein cyclin B, the regulatory phosphatase CDC 25 and the inhibitor CDC 6 known to be involved in the regulation of active CDK 1 and cyclin B complexes. The experimental data lead us to formulate a new hypothesis that CDC 6 slows down the activation of inactive complexes of CDK 1 and cyclin B upon mitotic entry. Our mathematical model, based on mass action kinetics, provides a possible explanation for the experimental data. We claim that the dynamics of active complexes CDK 1 and cyclin B have a similar nature to diauxic dynamics introduced by Monod in 1949. In mathematical terms we state it as the existence of more than one inflection point of the curve defining the dynamics of the complexes.

1. Introduction

The mitotic cell cycle is an ordered sequence of events, grouped into four phases: , , and , during which the eukaryotic cell doubles its content, and divides into two daughter cells. The classical models for studies of cell cycle molecular machinery are oocytes and early embryos. These have the distinguishing property of being transcriptionally silent. This implies that the molecular machinery governing oocyte maturation and early embryo development is based on the maternal information accumulated during oocyte growth. While many different proteins regulate the progression of the cell cycle and the transitions between cell cycle phases, the major regulatory mechanisms are based on similar processes in all phases. Two main classes of proteins involved in cell cycle control are cyclins and enzymes called cyclin dependent kinases—. During individual phases a specific cyclin accumulates in the cell, associates with an appropriate kinase and with the help of other enzymes activates the kinase/cyclin complex. The appropriate level of an active complex triggers the transition to the next phase of the cell cycle.

A protein complex responsible for the transition from to is formed from kinase and cyclin B, the latter being abbreviated as [3,4]. The complex of and , denoted by , can be in one of two states, inactive or active, abbreviated as and , respectively [3]. is dephosphorylated by active phosphatase, denoted by (inactive phosphatase is denoted by ). has the ability to pull away the phosphoryl group from two amino acids Tyr15 and Thr14 of the complex, making it active. This dephosphorylation results in the creation of . induces a cascade of phosphorylation of numerous substrates that change the character of cellular proteins from interphase to mitotic. These changes—necessary for mitotic progression—modify structures such as the cytoskeleton, membranes and DNA (condensation). Moreover, activation of phosphatase occurs due to its interaction with active complexes resulting in very powerful positive feedback between and that governs the activation upon the entry into -phase.

Summarising, a subtle equilibrium between and enzymes is maintained at the beginning of the to -phase transition process. The association of molecules with constantly synthesised results in the formation and accumulation of . This complex is inactive because it is phosphorylated on Tyr15/Thr14 since remains inactive. It was believed for a long time that, at some point, a spontaneous activation of the first molecules of triggers a positive feedback between and . Active molecules of start activation. Active , denoted as , activates new molecules of . This in turn triggers a dramatic acceleration of biochemical events leading to full activation of the whole pool of complexes present in the cell. Moreover, recently Vigneron et al. [5] have shown that another complex containing , namely, , triggers the activation of .

A number of mathematical models have been proposed to describe and understand mitotic cell cycle progression. For instance, some have used the stochastic approach to capture non-deterministic aspects of the process and emphasised that noise is an important factor influencing cell dynamics [6,7,8,9,10,11,12]. The common approach, however, is the deterministic one that usually is based on systems of ordinary differential equations. Some authors considered large systems of that type in order to investigate transitions between cell cycle phases. Such systems describing the dynamics of different proteins and enzymes were usually analysed only numerically, e.g., [13,14,15,16,17,18,19,20]. Some models have considered the activation of , e.g., [21,22,23]. There are cell cycle models based not only on standard ordinary differential equations, but also on delay differential equations [24,25]. For example, Busenberg et al. [24] contains a rigorous analysis of the model. Interesting analytic results can be found in other papers [26,27,28,29,30]. We emphasise the analytic results obtained by Ferrell et al. who described the cell cycle using ordinary differential equations and proved the existence of oscillatory dynamics of the cell cycle, switch-like behaviour of activity of and a bistability system [31,32].

The purpose of our work is to deepen the understanding of the cell cycle process. We are particularly focused on the to -phase transition process and we aim at investigating the role of the protein in entering into the -phase of the cell cycle. , up to now, was known as an essential ATPase active in the phase and responsible for the initiation of DNA replication.

Recent experiments made in Xenopus laevis one-cell embryo cell-free extract suggest that has an important role in the delay of to -phase transition [1]. This experimental system, however, simplified, as no nuclear DNA is present and all genetic information, both proteins and RNA, are of maternal origin and were accumulated in oocytes before the embryo development was triggered, allowed us to analyse the core molecular mechanisms of the cell cycle. We formulate a new hypothesis that explains this delay in terms of diauxic-like activation of . In a different context, diauxic behaviour was introduced by Monod [2] in 1949. We propose a new mathematical model that captures this new hypothesis. We present the analysis and numerical simulations of this new model, suggesting how regulates the dynamics of activation upon -phase entry.

The structure of the paper is as follows. In Section 2 we present the material and methods used to obtain experimental results that inspired our subsequent research. In Section 3, at first, we introduce the biochemical model describing the basic events occurring when a cell enters mitosis. Then, we present the experimental results on the protein that motivate our work. We show both the data reprinted from El Dika et al. [1] in Figure 1 and original results in Figure 2. We formulate a new hypothesis that captures the role of in the process and the mathematical model corresponding to the biochemical one. Next, we present the numerical simulations of the proposed model and finally, Section 4 provides conclusions and directions for further research. In Appendix A we present the mathematical analysis of the presented model.

Figure 1

activity in the control extract containing physiological amounts of (a) and in the extract immunodepleted of (b). Note a slow and diauxic growth of activity in the control extract (a) and the very rapid activation in the absence of (b). Curves reprinted from El Dika et al. [1].

Figure 2

Differences in dynamics of activation curves in control extracts containing physiological amounts of . Two extreme examples are shown: Rapid activation taking 16 min (a) and slow activation taking 28 min (b). Note that the inflection points of the curves appear at different moments in relation to the maximum activity.

2. Material and Methods

2.1. Egg Collection and Activation

Spawn Xenopus laevis eggs were dejellied with 2% l-cysteine pH 7.81 in XB buffer (100 mM KCl, 1 mM MgCl2, 50 mM CaCl2, 10 mM HEPES and 50 mM sucrose pH 7.6). Next, they were washed in XB buffer, activated with 0.5 mg/mL calcium ionophore A23187 and extensively washed in XB.

2.2. Cell Free Extracts

Cytoplasmic extracts from calcium ionophore-activated one-cell embryos before the first embryonic mitosis were prepared according to El Dika et al. [1]. In short, embryos were cultured at 21 °C in XB buffer for 60–70 min postactivation, transferred into 5 mL ultraclearTM centrifuge tubes (Beckman Coulter, Roissy, France) in 0.5 mL of XB buffer containing 0.1 mM AEBSF, a protease inhibitor, at 4°. They were subjected to three consecutive centrifugations: The first short spin to remove XB excess and pack the embryos, the second 10,000× g spin at 4 °C for 10 min to separate the cell-free fractions, and the final 10,000× g clarification spin of the supernatant at 4 °C for 10 min. The supernatant was then incubated at 21 °C. Aliquots were taken out every 4 min and stored at −80 °C.

2.3. 1 Activity Measurements

Samples of cell-free extracts were diluted in MPF buffer supplemented with: 0.5 mM sodium orthovanadate, 5 μg/μL of leupeptin, aprotinin, pepstatin and chymostatin, 0.4 mg/mL H1 histone (type III-S), 1Ci [32P] ATP (specific activity: 3000 Ci/mmol; Amersham Biosciences, UK) and 0.8 mM ATP. After incubation at 30 °C for 30 min, phosphorylation reactions were stopped by adding Laemmli sample buffer and heated at 85 °C for 5 min. Histone H1 was separated by SDS-PAGE and incorporated radioactivity was measured by autoradiography of the gel using a STORM phosphorimager (Amersham Biosciences, Buckinghamshire, UK) followed by data analysis with ImageQuant 5.2 software.

2.4. 6 Immunodepletion

Immunodepletion of from egg extracts was carried out using AffiPrep Protein A beads (Sigma, USA) conjugated with the anti- or with the preimmune serum overnight in 4 °C; 200 mL of beads were washed four times with XB buffer (pH 7.6) and incubated with 400 mL of extracts. After 30 min of incubation at 4 °C, extracts were centrifuged, beads were removed and supernatant was recovered. Two consecutive runs of immunodepletion were required to remove 90% of , as shown in El Dika et al. [1].

3. Results

3.1. Biochemical Model and the New Hypothesis

concentration gradually increases during the phase (cf. Equation (5)). pairs with protein kinase creates an inactive (phosphorylated) complex— (cf. Equation (1)). Inactive complex upon its interaction with phosphatase becomes activated, thus the concentration of active complexes increases (cf. Equation (2)). Conversely, complex activates phosphatase causing the appearance of more (cf. Equation (3)). Summarising, and form a positive feedback loop. The phase begins when the concentration of active exceeds the threshold value.

Recent experimental studies provoke intriguing questions about the role of the protein in slowing down the activation of complexes. Figure 1 shows the concentration of (from a biochemical point of view, simply the activity) obtained on the basis of molecular experiments in two cases: (a) With ; and (b) without (after removal of from the experimental system) [1]. In the experimental setting with , one can notice a slower increase in the concentration of . Therefore, in the experimental system with the entry into mitosis is delayed. Our main goal is to explain the role of in the observed phenomenon.

The diauxic growth of activity was clearly noticed in previous studies of Xenopus laevis one-cell embryo cell-free extracts ([33]: Figure 1A bottom, Figure 2A right, Figure 3A right, [34]: Figure 2A bottom and [28]: Figure 1V). Furthermore, in our own research, we always observed the same type of behaviour of ([35]: Figures 1A, 2A, 3A and 6A, [36]: Figures 2A, 3A, 6A and 7A,B and [37]: Figures 6A and 7B). Moreover, the diauxic growth of activity is not an artefact due to the cell-free system because it was also observed in individual Xenopus laevis one-cell embryos ([38]: Figure 1A); however, it is more clear in the vegetal hemisphere where the activation is delayed and proceeds with lower dynamics ([38]: Figure 1B). The precise dynamics of the diauxic growth of , i.e., inflection times and slope of the curve, varies form one experiment to another. For this reason, the average curves showing the dynamics of activation upon the -phase entry in Xenopus laevis one-cell embryos do not preserve the diauxic character ([35]: Figure 4B), where the average curve of 16 independent experiments does not show any inflection points. In Figure 2 in the current paper we show two examples of the fast and slow growth of activity in two independent experiments illustrating this problem well. The average curve of these two experiments also does not show the inflection points clearly visible in each experimental curve (data not shown).

We hypothesise that binds to and creates a new complex, preventing from being activated by phosphatase (cf. Equation (4)). The resulting complexes constantly break down into and that constantly associate again. The more accumulate in the cell, the more complexes are activated by residual . The formation of prolongs a very slow increase in the appearance of complexes. A slowdown in increase is visible as the flattening of the curve in Figure 1a. The experimental data leads to a new hypothesis on the mutual interaction between and , which determines the dynamics of upon mitotic entry. Our mathematical model, based on the law of mass action, bolsters this hypothesis. We suggest that the dynamics of are similar to diauxic dynamics introduced by Monod [2]. In mathematical terms, we state it as the existence of more than one inflection point on the curve defining the dynamics of the complexes, cf. Figure 3. Indeed, in the present model, we observe three or four inflexion points.

Figure 3

The smoothed curve obtained on the basis of experimental data presented in Figure 1a. Red circles indicate approximate location of inflection points for the setting with upon -phase entry.

The second part of our hypothesis is that the reaction speed of and binding depends on active in a switch like mode. This means that when the concentration of is less than the concentration value, then the reaction speed of formation is low. When the is higher than the threshold value the reaction speed becomes much faster resulting in a two-step activation visible in biological experiments as an inflection of the activation curve of . Experimental data show clearly that this activation depends on the presence of [1].

Summarising, we consider the biochemical model that takes into account eight species, the descriptions of which are provided in Table 1, whereas the scheme of their mutual interactions is provided in Figure 4.

Table 1

Species of the biochemical model involved in Equations (1)–(5).

Species	Description
CDK1	cyclin-dependent kinase 1
CYCB	cyclin B
CDK1/CYCB_A	active complex of CDK and CYCB
CDK1/CYCB_N	inactive complex of CDK and CYCB
CDC25_A	active phosphatase CDC25
CDC25_N	inactive phosphatase CDC25
CDC6	cell division cycle 6 ATPase
CDK1/CYCB/CDC6	complex of CDK1/CYCB_N and CDC6

Figure 4

The schematic diagram of the considered system. Colours of arrows and dots correspond to colours of Equations (1)–(5). For simplicity we do not consider the potential marginal separation of the complex into and . On the diagram we indicate this by “*”.

We consider the following five reactions (i.e., Equations (1)–(5)).

We want to emphasise that Equations (1)–(3) and (5) correspond to the current state of knowledge. Equation (4) reflects the new hypothesis and is our contribution to understanding the phenomenon. In summary, taking Equation (4) into account is the first part of the new hypothesis. The speed of Equation (4) described by the function f is the second part of the hypothesis.

3.2. Mathematical Model

Assuming mass action kinetics for Equations (1)–(5) we transform the biochemical model into the system of eight ordinary differential equations (ODEs) with the following notation

We have where are positive parameters and . The function f is the Hill function that describes switch-like behaviour, where is a positive coefficient, k is a Hill coefficient, is the threshold value of the switch and is the basic rate when . The function f describes the reaction rate of associated with resulting in the formation of complexes (cf. Equation (4)). In the system of ordinary differential equations, Equation (6) appears in the 7th equation describing the dynamics of and, due to the law of mass action, in the 3rd and 6th equations describing the dynamics of and , respectively. The process of formation seems to be highly nonlinear and we assume its rate to be dependent. There exists a similar mechanism governing interactions between and in -phase. If is low then the majority of is not phosphorylated. However, with an increase of more phosphorylated appears in the cell [39]. The function f is bounded as it plays the role of a rate coefficient. The typical way of modelling such a nonlinear dependence is based on the Hill function, see, e.g., [40].

Taking into consideration the biological constraints, we propose the following initial data

Equation (6) have the following conservation laws where denote constants given at the initial time. In Appendix A we provide the mathematical analysis of the model.

We provide the standard non-dimensionalisation of Equation (6). In other words we relate all considered variables to their characteristic values. With the substitution and omitting the stars for simplicity, we obtain

By Equation (8) it follows

By Equation (7) we obtain

From the mathematical analysis presented in Appendix A we deduce that if the system contains even a small amount of or then and converge to full activation. This result is consistent with biological observations, because if the initial concentration of or is positive then the positive feedback loop starts and the biological system tends to its equilibrium state (called ) defined by the maximal concentrations of and . If the initial concentrations of or are equal to zero, then the positive feedback loop does not start and the biological system tends to another equilibrium state (called ) defined by the concentrations of and equal to 0. Small perturbations of the initial concentrations from zero to positive values change the equilibrium points, and this is the biological reason for being unstable and being asymptotically stable.

We note that a further simplification of the reduced model, Equation (A3), considered in Appendix A is reasonable. For example taking and we may reduce this system to a system of three equations

3.3. Numerical Simulations

To carry out the numerical simulations we use the Runge-Kutta 4th order method provided by Matlab. Parameters values used to carry out the numerical simulations are given in Table 2. Figure 5 shows the concentrations of , , , . The most interesting curve is , where we observe three inflection points. The concentrations of species containing are shown in Figure 6 with the concentration of set to zero. Figure 7 shows the difference in activation: The timing and dynamics of the activation of in the presence and absence of . When is present the activation has more than one inflection point, and mitosis starts later, whereas when is absent, the activation is fast. In Figure 5 and Figure 7 we observe diauxic-type behaviour for the curve of . According to our hypothesis, this is related to the mutual interaction between and . We link this kind of behaviour with the existence of multiple (three or four in this case) inflection points in the curve of . The rigorous investigation of this fact leads to the analysis of behaviour of the second derivative of and more precisely its number of zeros. Figure 8 presents the graphs of the second derivative of obtained for the reduced system of three equations, Equation (13).

Table 2

Values of parameters used in simulations.

Parameter	Value	Parameter	Value	Parameter	Value
α_1	1	δ	4	ε_x_a	0.001
α_2	30	k	20	ε_x_n	0.001
α_3	1	ν	8	ε_y_a	0
α_4	7	ν_th	0.25	ε_w	0
γ	1	ω	0.6

Figure 5

Concentration of , , in the presence of .

Figure 6

Concentration of , , in the absence of .

Figure 7

Comparison between concentration of in the absence and presence of . Solid line—system with ; dotted line—system without .

Figure 8

Graphs presenting second derivatives of showing the number of zeros, which indicates the number of inflection points. (a) corresponds to the case with the second derivative starting from a positive value and having three zeros. (b) corresponds to the case with the second derivative starting from a negative value and having four zeros.

We may note that the numerical result given in Figure 8 has a rigorous nature as is visible by a careful estimation of the error of Matlab approximation. For example, considering Equation (13), we may provide (following the idea of [41]) the detailed analysis of the error where is the value of the true solution at point and is the approximation of the solution at point , showing that

Taking the error is smaller than the variation in the second derivative. Moreover, the rounding error can be neglected because the machine epsilon (see [42]) is sufficiently small compared with h. The details of the estimation are not reported here. This leads to the conclusion that the result stating the number of zeros of the second derivative has a rigorous nature and there the number of inflection points of the variable is either three or four, which give diauxic behaviour.

4. Discussion

The proposed model captures the most important characteristics of the diauxic growth of activation observed in biochemical experiments. Based on our previous experimental results [1] we claim that is the most important factor which causes the inflection of the activation curve. We have shown for the first time that is an inhibitory protein acting on during -phase. Making use of our modelling setting, we hypothesise that binds to forming . formation results in slower activation of and consequently a delayed entry into -phase. From a biochemical perspective our results, both experimental and modelling, are particularly interesting because the inhibitory effect of on activation during -phase was not shown previously.

The second part of our hypothesis stands for the switch-like dependence of the reaction rate of binding to resulting in the formation of (Equation (4)). Our assumption that the mentioned reaction rate depends on provides a good qualitative explanation of the observed diauxic dynamics of activation. Further biological research is needed to investigate what molecular modification is necessary for this switch-like pattern of activation. We can postulate that , at some threshold, phosphorylates triggering the abrupt increase in affinity to . From a more general perspective, the slow rate of activation is very likely important for the physiological course of mitotic processes such as chromatin condensation or spindle formation in such a large cell as the Xenopus laevis one-cell embryo.

One of the main goals of the paper was to describe mathematically the diauxic behaviour of activation in the presence of protein. However, this kind of approach has a wide spectrum of use and may be applied to a large variety of problems. Usually such complex biological systems are very difficult to treat rigorously from a mathematical point of view. The analysis of the error gives a chance for a rigorous statement based on the numerical simulations. We leave the details of this approach for a forthcoming paper.

Our results may affect the understanding of the process of cancerogenesis since and its interactions with play an important role in mitotic regulation and in cancer etiology [43,44]. The role in -phase regulation is not limited only to the mitotic cell cycle as shown in the current paper and in El Dika et al. [1], but also to the meiotic regulation in oocytes [45,46,47,48]. The requirement of for the meiotic spindle formation in mice and Xenopus laevis oocytes suggests that it can also be involved in mitotic spindle formation. The diauxic growth of activity determined by may be in relation to the proper dynamics of spindle assembly not only through the fine tuning of microtubule dynamics, but also by the proper coordination with other players like actin filaments [49,50].