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 CDC6 protein in the entry of cells into mitosis. This has not been considered in the literature so far. Recent experiments suggest that CDC6, upon entry into mitosis, inhibits the appearance of active CDK1 and cyclin B complexes. This paper proposes a mathematical model which incorporates the dynamics of kinase CDK1, its regulatory protein cyclin B, the regulatory phosphatase CDC25 and the inhibitor CDC6 known to be involved in the regulation of active CDK1 and cyclin B complexes. The experimental data lead us to formulate a new hypothesis that CDC6 slows down the activation of inactive complexes of CDK1 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 CDK1 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.


Introduction
The mitotic cell cycle is an ordered sequence of events, grouped into four phases: G 1 , S, G 2 and M, 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-CDKs. 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. biochemical model describing the basic events occurring when a cell enters mitosis. Then, we present the experimental results on the CDC6 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 CDC6 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.

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 MgCl 2 , 50 mM CaCl 2 , 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.

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 ultraclear TM 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.

CDK1 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), 1µCi [γ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.

CDC6 Immunodepletion
Immunodepletion of CDC6 from egg extracts was carried out using AffiPrep Protein A beads (Sigma, USA) conjugated with the anti-CDC6 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 CDC6, as shown in El Dika et al. [1].

Biochemical Model and the New Hypothesis
CYC B concentration gradually increases during the G 2 phase (cf. Equation (5)). CYC B pairs with protein kinase CDK1 creates an inactive (phosphorylated) complex-CDK1/CYC B N (cf. Equation (1)). Inactive complex CDK1/CYC B N upon its interaction with phosphatase CDC25 A becomes activated, thus the concentration of active complexes CDK1/CYC B A increases (cf. Equation (2)). Conversely, complex CDK1/CYC B A activates phosphatase CDC25 N causing the appearance of more CDC25 A (cf. Equation (3)). Summarising, CDC25 A and CDK1/CYC B A form a positive feedback loop. The M phase begins when the concentration of active CDK1/CYC B A exceeds the threshold value.
Recent experimental studies provoke intriguing questions about the role of the CDC6 protein in slowing down the activation of CDK1/CYC B N complexes. Figure 1 shows the concentration of CDK1/CYC B A (from a biochemical point of view, simply the CDK1 activity) obtained on the basis of molecular experiments in two cases: (a) With CDC6; and (b) without CDC6 (after removal of CDC6 from the experimental system) [1]. In the experimental setting with CDC6, one can notice a slower increase in the concentration of CDK1/CYC B A . Therefore, in the experimental system with CDC6 the entry into mitosis is delayed. Our main goal is to explain the role of CDC6 in the observed phenomenon. The diauxic growth of CDK1 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 CDK1 ( [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 CDK1 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 CDK1 activation is delayed and proceeds with lower dynamics ( [38]: Figure 1B). The precise dynamics of the diauxic growth of CDK1, 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 CDK1 activation upon the M-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 CDK1 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 CDC6 binds to CDK1/CYC B N and creates a new CDK1/CYC B/CDC6 complex, preventing CDK1/CYC B N from being activated by CDC25 A phosphatase (cf. Equation (4)). The resulting CDK1/CYC B/CDC6 complexes constantly break down into CDC6 and CDK1/CYC B N that constantly associate again. The more CDK1/CYC B N accumulate in the cell, the more CDK1/CYC B N complexes are activated by residual CDC25 A . The formation of CDK1/CYC B/CDC6 prolongs a very slow increase in the appearance of CDK1/CYC B A complexes. A slowdown in CDK1/CYC B A 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 CDC6 and CDK1/CYC B N , which determines the dynamics of CDK1/CYC B A upon mitotic entry. Our mathematical model, based on the law of mass action, bolsters this hypothesis. We suggest that the dynamics of CDK1/CYC B A 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.
The second part of our hypothesis is that the reaction speed of CDK1/CYC B N and CDC6 binding depends on active CDK1/CYC B A in a switch like mode. This means that when the concentration of CDK1/CYC B A is less than the concentration value, then the reaction speed of CDK1/CYC B/CDC6 formation is low. When the CDK1/CYC B A is higher than the threshold value the reaction speed becomes much faster resulting in a two-step CDK1 activation visible in biological experiments as an inflection of the activation curve of CDK1. Experimental data show clearly that this activation depends on the presence of CDC6 [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.  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.  (1)-(5). For simplicity we do not consider the potential marginal separation of the complex CDK1/CYC B N into CDK1 and CYC B. On the diagram we indicate this by " * ".

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 where α 1 , α 2 , α 3 , α 4 , β, δ are positive parameters and f (x) = ω + νx k v k th +x k . The function f is the Hill function that describes switch-like behaviour, where ν is a positive coefficient, k is a Hill coefficient, v th is the threshold value of the switch and ω is the basic rate when x a = 0. The function f describes the reaction rate of CDK1/CYC B N associated with CDC6 resulting in the formation of CDK1/CYC B/CDC6 complexes (cf. Equation (4)). In the system of ordinary differential equations, Equation (6) appears in the 7th equation describing the dynamics of CDK1/CYC B/CDC6 and, due to the law of mass action, in the 3rd and 6th equations describing the dynamics of CDK1/CYC B N and CDC6, respectively. The process of CDK1/CYC B/CDC6 formation seems to be highly nonlinear and we assume its rate to be CDK1/CYC B A dependent. There exists a similar mechanism governing interactions between CDK1 and CDC6 in S-phase. If CDK1/CYC B A is low then the majority of CDC6 is not phosphorylated. However, with an increase of CDK1/CYC B A more phosphorylated CDC6 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 K CDK1 , K CDC25 , K CDC6 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 obtaiṅ x = −α 1 xc , x a = α 2 x n y a , x n = α 1 xc − α 2 x n y a − α 4 f (x a )x n z + δw , y a = α 3 x a y n , y n = −α 3 x a y n , By Equation (8) it follows x a + x n + x + w = 1, y a + y n = 1, z + w = γ .
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 CDK1/CYC B A or CDC25 A then CDK1/CYC B A and CDC25 A converge to full activation. This result is consistent with biological observations, because if the initial concentration of CDK1/CYC B A or CDC25 A is positive then the positive feedback loop starts and the biological system tends to its equilibrium state (called S 2 ) defined by the maximal concentrations of CDK1/CYC B A and CDC25 A . If the initial concentrations of CDK1/CYC B A or CDC25 A are equal to zero, then the positive feedback loop does not start and the biological system tends to another equilibrium state (called S 1 ) defined by the concentrations of CDK1/CYC B A and CDC25 A 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 S 1 being unstable and S 2 being asymptotically stable.
We note that a further simplification of the reduced model, Equation (A3), considered in Appendix A is reasonable. For example taking x = 0 and c = 1 we may reduce this system to a system of three equationsẋ

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 CDK1, CDK1/CYC B A , CDK1/CYC B N , CDK1/CYC B/CDC6. The most interesting curve is CDK1/CYC B A , where we observe three inflection points. The concentrations of species containing CDK1 are shown in Figure 6 with the concentration of CDC6 set to zero. Figure 7 shows the difference in activation: The timing and dynamics of the activation of CDK1/CYC B A in the presence and absence of CDC6. When CDC6 is present the activation has more than one inflection point, and mitosis starts later, whereas when CDC6 is absent, the activation is fast. In Figures 5 and 7 we observe diauxic-type behaviour for the curve of CDK1/CYC B A . According to our hypothesis, this is related to the mutual interaction between CDC6 and CDK1/CYC B N . We link this kind of behaviour with the existence of multiple (three or four in this case) inflection points in the curve of CDK1/CYC B A . The rigorous investigation of this fact leads to the analysis of behaviour of the second derivative of CDK1/CYC B A and more precisely its number of zeros. Figure 8 presents the graphs of the second derivative of x a obtained for the reduced system of three equations, Equation (13).

Parameter Value Parameter Value Parameter Value
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 y(t n ) is the value of the true solution at point t n and y n is the approximation of the solution at point t n , showing that |E n | < 4.5 · 10 10 · h 4 .
Taking h = 1 1000 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 x a is either three or four, which give diauxic behaviour.

Discussion
The proposed model captures the most important characteristics of the diauxic growth of CDK1 activation observed in biochemical experiments. Based on our previous experimental results [1] we claim that CDC6 is the most important factor which causes the inflection of the CDK1 activation curve. We have shown for the first time that CDC6 is an inhibitory protein acting on CDK1 during M-phase. Making use of our modelling setting, we hypothesise that CDC6 binds to CDK1/CYC B N forming CDK1/CYC B/CDC6. CDK1/CYC B/CDC6 formation results in slower activation of CDK1 and consequently a delayed entry into M-phase. From a biochemical perspective our results, both experimental and modelling, are particularly interesting because the inhibitory effect of CDC6 on CDK1 activation during M-phase was not shown previously.
The second part of our hypothesis stands for the switch-like dependence of the reaction rate of CDK1/CYC B N binding to CDC6 resulting in the formation of CDK1/CYC B/CDC6 (Equation (4)).
Our assumption that the mentioned reaction rate depends on CDK1/CYC B A provides a good qualitative explanation of the observed diauxic dynamics of CDK1 activation. Further biological research is needed to investigate what molecular modification is necessary for this switch-like pattern of CDK1 activation. We can postulate that CDK1, at some threshold, phosphorylates CDC6 triggering the abrupt increase in CDC6 affinity to CDK1/CYC B N . From a more general perspective, the slow rate of CDK1 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 CDK1 activation in the presence of CDC6 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 CDC6 and its interactions with CDK1 play an important role in mitotic regulation and in cancer etiology [43,44]. The CDC6 role in M-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 CDC6 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 CDK1 activity determined by CDC6 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].

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Mathematical Analysis of the Model
The global existence, uniqueness, positiveness and boundedness of the solutions follow directly from the form of Equation (10). Keeping in mind Equation (12), we obtain for all t > 0. From Equation (11) we obtain Using Equation (A2) we reduce Equation (10) tȯ Referring to equilibrium points we formulate the following proposition.
Referring to the stability of the equilibrium points.
Proposition A2. The equilibrium points S 1 and S 2 are unstable and asymptotically stable, respectively.
Next we study the stability of the point S 2 . The corresponding Jacobi matrix is J(S 2 ) is a block matrix; thus, we have where R 1 = − α 1 , The characteristic polynomial w R 2 of matrix R 2 is: Each coefficient of polynomial w R 2 is negative, so each root has a negative real part. Matrices R 1 and R 3 have eigenvalues −α 1 and −1, respectively, which are negative. Each eigenvalue of matrix J(S 2 ) has a negative real part, so the stationary state S 2 is asymptotically stable.