A quantitative model of the initiation of DNA replication in Saccharomyces cerevisiae predicts the effects of system perturbations

Description of model components

The model describes sixteen molecular species (twelve of which are dynamically independent) and depends on twenty-four parameters, which characterize the rates of seventeen biochemical processes (protein expression and degradation, complex association/dissociation, and transport across the nuclear membrane). The model describes the following molecular species (Figure 1).

The MCM species corresponds to dimers of Mcm2-7 heterohexamers, as two complexes are loaded at each origin. Similarly, the CDC45 species corresponds to a dimer, as described by Bowers et al. [11]. We describe concentrations in units of molecules per cell.

Reaction events

In choosing reaction kinetics, we balanced the complexity of the model against its ability to adequately describe the behaviour of the overall system. We limited our description of initiation to the interactions between the pre-RC and replisome proteins that we found to be the essential core of the network (e.g. Dbf4 representing the Dbf4-Cdc7 complex, discussed below). As a result, certain processes were combined into single events, some reactions were presumed irreversible, and only some reaction rates were presumed to have non-linear kinetics.

Except for RC7, phosphorylation states are not explicitly described, as we have no data for the individual phosphorylation events. This is acceptable for our purposes as the lumped function of CDK in each case is consistent with a scenario where the effect of CDK is proportional to its concentration (i.e. [Clb5]). Additionally, processes that involve multi-protein complexes are represented by a single member – one CMG (Cdc45·Mcm2-7·GINS) complex stabilizes DNA polymerase at each replication fork. Of the three protein factors it is comprised of, Cdc45 is limiting. Although MCM is also included in the GINS complex, we model both MCM and Cdc45 as separate species. Dbf4 represents the Dbf4-Cdc7 kinase complex and Mcm2 represents the Mcm2-7 helicase. Although the protein factors Cdc45, Dbp11, Sld2, Sld3 and GINS interact to facilitate formation of the pre-initiation complex at origins, we model only Cdc45, which is the limiting factor in the CMG complex [55]; ultimately the number of forks fired (described by our model) is dependent on the Cdc45 concentration. We take Dbf4, which is the limiting regulatory subunit of Cdc7, as representative of active DDK, which is one of the limiting factors in replication initiation [56]. Mcm2 is used to represent MCM complexes; the Mcm2 concentration has been reported to approximate the number of total complexes per cell in an asynchronous population [57, 58]. The replication complexes in our models exist only on chromatin and therefore represent the activity of these proteins at the DNA as opposed to soluble complexes.

The network shown in Figure 1 includes both reversible and irreversible reactions as indicated. Association/dissociation reactions are considered reversible, in accordance with a dynamic pre-RC/pre-IC loading mechanism as described above. In most cases, phosphorylation events are modeled as irreversible, in the absence of identified countervailing enzymes. We found that it was sufficient to describe most reaction rates by mass action kinetics. In cases where saturation occurs (the nuclear import of MCM·Cdt1, v₆, and the association of Cdc45 with ORC, v₁₃), we employed Michaelis-Menten kinetics. To simplify the description of the phosphorylation of ORC by CDK (RC7), we do not describe phosporylation and dephosphorylation explicitly, but combine them into a single dephosphorylation event whose rate is inversely proportional to the level of CDK (v₁₇). We introduce cooperativity in this mechanism to account for multiple phosphorylation events [32] or an additional inhibitory CDK-Orc6 binding mechanism [34, 35].

The establishment of replication complexes in our model reflects the sequential binding of proteins that constitute the pre-replicative complex. In some cases the association and dissociation of pre-RC components is reversible. We treat the loading and maintenance of Mcm2-7 helicase complexes as a dynamic process, which is dependent on the concentrations of the factors ORC, Cdc6, Cdt1 and Mcm2-7 itself. A mechanistic model for the dynamic assembly of pre-RCs was first described by the Bell lab [59]. The requirement of pre-RC factors for maintenance of helicase-loaded origins in late G1 has been further demonstrated by work from these researchers as well as our group [8, 9, 60].

Network and differential equations

where RC_Total, CDT1_Total, MCM_Total, and CDC45_Total are the fixed total number of origins, Cdt1 molecules, Mcm2-7 complexes, and Cdc45 dimers, respectively. These four factors have been shown to be present at constant levels throughout the cell cycle [7, 61–64]. The value for RC_Total used in our model is 332, as described in [65].

System inputs

The biological network responsible for the initiation of DNA replication does not oscillate autonomously; it displays periodic behaviour when driven by periodic signals from the cell cycle. Likewise, our model displays oscillations only when driven by periodic forcing input. In order to facilitate the combination of our model with the cell cycle model of Chen et al. [45], we used the simulated profiles of Clb5 and Cdc20 from their model as periodic inputs to ours. Cdc20 mediates the degradation of Dbf4 (reaction v₄). Clb5 is responsible for Cdc6 degradation (v₂), loading of Cdc45 (v₁₃), nuclear export of free MCM (v₁₆) and Cdt1 (v₈), and phosphorylation of Orc2 and Orc6 (v₁₇). We converted the time-varying profiles of Clb5 from the Chen et al. [45] model to molecules-per-cell units using the genome-wide GFP tagging and localization experiments described in [66, 67]. The profile of Cdc20 was similarly obtained by scaling to cellular abundance levels reported in another study – while Cdc20 has been determined to peak at 2200 copies in a haploid cell, the functional APC^Cdc20 level can be estimated by considering the APC cyclosome subunit Cdc27 [68, 69]. This was reported in different studies to be 593 mol/cell in an asynchronous population [67] and at its maximal value of 750 mol/cell in metaphase [68].

Data acquisition

Data for Cdc45 and Cdc6 levels were obtained from individual isogenic strains in which the open reading frame of the corresponding gene was fused to a sequence encoding a 13Myc epitope tag [70]. In Figure 2, a representative western blot for Cdc45-Myc is shown (panel A), with the corresponding FACS analysis (panel B). The levels of Mcm2 were determined using an anti-Mcm2 antibody. In each of our time course experiments, cells were first arrested in late G1 phase with the mating pheromone α-factor and then released synchronously into the cell cycle, as described in Methods. From the literature, we used time course data for chromatin-bound and soluble Dbf4 and Mcm2 (to supplement our in vivo data) from Pasero et al. [71] and quantitation of the nuclear fraction of Cdt1 from Tanaka and Diffley [7]. In order to convert relative measures of protein abundance to molecule-per-cell numbers, we used scaling factors obtained from the database provided by Ghaemmaghami et al. [67]. The data is shown along with a best-fit simulation in Figure 3. Raw timecourse data can be found in Additional files 12 and 3. Our calculation of molecule-per-cell estimates is demonstrated in Additional file 4. The fits in Figure 3 represent our best solution to a trade-off between quality of fit and model complexity. We explored the effect of adding additional species and parameters. These additional features could, in some cases, provide minor improvements to the fit, but our confidence in the parameter estimates (discussed in the following section) suffered as the complexity of the model grew.

While our model reports proteins and complex concentrations in absolute units of molecules/cell, our accuracy regarding these values is limited to by the literature-reported cellular abundances for the various protein factors. We have used discretion when inconsistencies arise, choosing the reported values that are best justified by multiple studies. The network dynamics are based on the relative protein abundances over the time-course (indeed, many modelling efforts describe concentration changes in arbitrary units). Consequently, changes in the global protein level for a particular factor do not affect the dynamics; instead, they impact our predictions of absolute molecule/cell counts and parameter values. This is an important consideration regarding the conversion of densitometry readings to absolute values, as the overall levels are ultimately determined by a literature-derived scaling factor. Future experiments will undoubtedly result in improved estimates of protein abundances. These can be easily incorporated into the model by scaling the parameter values accordingly (with no direct impact on system dynamics).

Parameter calibration

Parameter values that were well constrained by the data include those specifying the rates of production, degradation and association of Cdc6 (k₁, k₂, k₅) and Dbf4 (k₃, k₄, k₁₂) as well as the rate of origin firing (k₁₄). This reflects the strong reliability of our data for these two protein factors as well as for the proteins that form the replication complex (RC6) that gives rise to active forks.

Parameters values in which we have low confidence include those that govern the loading of MCM by Cdt1 (v₇), Cdc6 dissociation from RC3 (k_5r), and the phosphorylation of ORC (k₁₇). The reversible dissociation of Cdc6 is needed to accurately fit the data and there is no evidence suggesting that ORC-Cdc6 binding is irreversible. Nevertheless, it is clear that experimental observations specific to this process are required to more precisely estimate this parameter value. The reaction whereby the MCM·Cdt1 species loads the MCM complex (v₇) is extremely transient [10]. Provided parameter k₇ is sufficiently large, the kinetics of this reaction will be rapid enough to fit the data. Consequently, the data cannot support a precise estimate of the parameter value. This observation suggests that MCM loading is an extremely rapid biochemical step in pre-RC assembly. It may point to a role for Cdt1 in repeatedly targeting MCM complexes to origins throughout G1. Such a phenomenon is consistent with the requirement for a dynamic loading mechanism that ensures pre-RC fidelity up until the G1/S transition. Finally, the phosphorylation of ORC (characterized by k₁₇) contributes to the prevention of repeated origin firing. However, this mechanism has not been well characterized, and our data is unable to accurately constrain the specifics of this process.

The kinetic rates in this network have not been the subject of prior experiments, but previous reports of protein half-lives are consistent with our predicted parameter values. Drury et al. [31] estimated that Cdc6 is reduced below the point of detection within 5 minutes of S-phase entry, corresponding to a half-life no longer than 1.5 min. Similarly, Cheng et al. [20] reported that Dbf4 is reduced below visible levels within 10 minutes by the APC-dependent pathway, indicating a half-life no longer than 3 min. Our model-based predictions of degradation rates correspond to half-lives of 1 min. and 2.5 min. for Cdc6 and Dbf4 respectively, in good agreement with these earlier findings.

Figure 4 shows the simulated model behaviour for the best-fit parameter set. Some replication complex species – RC2, RC4 and RC5 – are extremely transient. Their low levels of abundance are shown separately from other RCs, on an appropriate scale. Simulations were carried out in Matlab (code available from the authors upon request).

Perturbations

Our initial explorations of the model revealed that the network’s behaviour is particularly sensitive to the abundance of Dbf4 and Cdc6 and relatively insensitive to the level of Cdt1. We investigated the effects of perturbations by simulating reductions in Dbf4, Cdt1 and Cdc6 (Figure 5) in our model. When the Cdc6 production rate (v₁) was reduced to 10% of its nominal (wild-type) value, persistence of the RC1 complex was observed. Similarly, when the Dbf4 production rate (v₃) is reduced by the same relative amount, an accumulation of RC3 occurs. In both cases, the perturbation interferes with pre-initiation complex assembly and blocks the system at the nearest previous persistent RC state; RC4 is not persistent since the unloading of MCM causes a rapid transition back to RC1. It is worth noting that because MCM can dissociate from ORC (v₂₂), RC4 represents a complex containing MCMs that will be functionally incorporated into replication forks as opposed to those that loosely associate with origins. Because the timing of our model is fixed, the various state concentrations (RC levels) indicate the progression from licensing to firing. A reduction in the FORK species compared to the wild-type case suggests a slow-down in S-phase because fewer origins are firing within the prescribed time. Using the peak abundance of the FORK species as a measure of replicative efficiency, we saw significant reductions in both simulated knock-downs (by 68% for Dbf4 and 73% for Cdc6, Figure 5B, C). Conversely when we simulated the reduction of Cdt1 abundance to 10% of nominal values, origin firing was only reduced by 23%, suggesting that the network is relatively refractory to depletion of Cdt1 (Figure 5D). Additional file 5 shows the levels of the various model components for each perturbation.

To investigate the accuracy of these mathematical predictions, we carried out corresponding wet lab depletion experiments. Reducing Dbf4 or Cdc6 concentrations in yeast cells to roughly 90% below normal endogenous levels resulted in a rapid G1 phase arrest, evident after 2 hours of depletion, as judged by FACS analysis indicating the accumulation of cells with 1 C (unreplicated) DNA content (Figure 6). In contrast, a corresponding depletion of Cdt1 had no appreciable effect, and DNA replication defects were still only minimally evident after 8 hours of further reduction. Thus, our in silico simulations using our nominal parameter set were predictive of in vivo perturbations. These experiments were used to validate the model; they were not used for calibration.

The insensitivity to perturbations in Cdt1 levels is consistent with its apparent excess relative to origins [67], although the number of Cdt1 molecules that act at each origin has not yet been characterized. Moreover, the mechanism by which Cdt1 aids in recruiting the helicase molecules to pre-RCs is extremely transient [10].

While many factors are limiting, the system appears to be highly sensitive to the levels of Cdc6. Due to its low abundance relative to MCM and Cdt1, even a moderate depletion of Cdc6 significantly alters the dynamics of pre-RC loading. The same is true for Dbf4, although in this case its role in activation of the Cdc7 kinase renders the system highly sensitive to its concentration; firing cannot occur without the Dbf4-Cdc7 complex. Since Dbf4 is, like Cdc6, limiting, flow through the network is blocked when the kinase does not reach a threshold level. Additionally at limiting levels, the number of replication forks produced by our model is significantly reduced, consistent with in vivo reports from the literature showing a lengthening of S-phase [15].

Further validation of our model comes from comparison with additional in vivo experiments reported in the literature. Jones et al. [73] showed that the interaction between the MCM complex and Dbf4 was reduced to half its wild-type level when a Dbf4 domain that binds Mcm2 was mutated, impairing S-phase progression. We mimicked this effect by reducing the rate of association of Dbf4 with RC4 (k₁₂) by 50%, leading to a similar result (compare Figure 7 panels A and B). Similarly, Zou et al. [74] reported that the cdc45-1 mutant shows an aberrant growth phenotype at the non-permissive temperature. This is thought to be due to a disruption of Cdc45’s ability to interact with MCM and ORC (RC6). As shown in Figure 7C, by reducing the rate of Cdc45 interaction with RC6 (k₁₃) by 50%, a marked reduction in the peak abundance of the FORK species results, indicative of a slower S-phase, as observed when the mutant was grown at the non-permissive temperature. The actual reduction in Cdc45’s association with the pre-RC due to conformational changes in the mutant might be even more pronounced than a 2-fold reduction. In any case, our simulation is consistent with Cdc45’s origin-initiation role being compromised by impairing its ability to interact with its ligands to form the CMG complex.

Linking our DNA initiation model to a previously established cell cycle model

Our model of the initiation of DNA replication only displays oscillatory behaviour when forced with periodic signals from the cell cycle. By choosing to incorporate signals that correspond to species in the cell cycle model of Chen et al. [45], we were able to merge the two models in a straightforward manner. In the Chen model, the initiation of DNA replication is represented by a single lumped state variable, called ORI. At the beginning of the cell cycle, ORI has value zero. Its rate of growth depends linearly on Clb5. When it reaches a threshold value, DNA synthesis is presumed to have begun, and triggers an increase in the value of the parameter k_mad2 (activity level of the Mad2 protein) leading to an inactivation of Cdc20, which is required for mitotic exit. This Mad2-dependent inhibition of Cdc20 represents the spindle assembly checkpoint [75], ensuring that cells with replicated DNA do not complete mitosis without properly aligning the chromosomes. When chromosomes have properly aligned on the metaphase plane k_mad2 drops and Cdc20 promoted exit from mitosis. In our model, the level of DNA synthesis is represented by the FORK species. To merge the two models, we removed the ORI state from the Chen model, and instead used the FORK species to trigger the change in k_mad2, as detailed in Methods.

Besides “closing the loop” between the two models by incorporating two-way inter-model signalling (involving Clb5, Cdc20, and FORK), we had to deal with a single shared species: both models describe the dynamics of Cdc6. We arrived at a merged description of Cdc6 behaviour by incorporating the dynamics of recognition complex association and dissociation into the Chen model’s formulation of Cdc6 behaviour (details are described in Methods). The resulting combined model behaves only marginally differently from either model in isolation, as shown in Additional file 6: Figure S6 and Additional file 7: Figure S7).

Replication metrics of Brümmer et al

The recently published model of Brümmer et al. [54] also describes the network responsible for the initiation of DNA replication. The 51 free parameters of that model were chosen by a combination of fitting and optimization. The authors used literature-derived data to fix 28 of the kinetic parameters. The remaining 23 free parameters were not fit to data, but were selected through a procedure that optimized the coherence of origin firing and minimized re-replication (selected as hypothetical goals of evolutionary ‘design’). While it is impossible to assess the accuracy of the parameter values obtained from this procedure, the resulting idealized model provides a useful starting point for examining how the network structure constrains the system behaviour.

The model of Brümmer et al. focuses on early origin firing and so represents the mechanics of firing at the start of S-phase. In contrast, our model describes firing dynamics throughout S-phase in order to fit into the broader context of the cell cycle [76–79]. The parameter set driving our system is not filtered to retain only those that produce replication dynamics consistent with coherent firing just at the G1/S boundary. Rather, the parameters are specified by the actual cellular concentrations of the active protein factors generating replication forks. While both models incorporate the important role of CDK, Brümmer et al. emphasize the multi-site phosphorylations of several factors involved in mechanisms that minimize potential re-replication. To this end, Brümmer et al. employed a metric to assess re-replication. Their idealized model exhibits 0.0028 re-replication events per cycle. Applying the same measure to our bottom-up model yields 0.36 re-replication events per cycle (although that can readily be reduced by modifying our parameters from their best-fit values). The near-zero value obtained by Brümmer et al. is close to their idealized target of zero. Both estimates are consistent with the belief that re-replication occurs in wildtype cells, but at an extremely low rate [79]. Because the nature of the dephosphylation of ORC (RC7→RC1 transition) remains uncharacterized, we use a conservative estimate of the number of ORC phosphorylation sites. Increasing this number by twofold, consistent with the number of CDK target residues on ORC [32] reduces re-replication to a value on the same order of magnitude as Brümmer’s value. Thus, both models effectively deal with representing control and prevention of rereplication.

Yeast strains and transformants

Cell cycle timecourse

Cells were grown in YPD medium (1% Bacto-yeast extract, 2% Bacto-peptone, 2% dextrose) to exponential phase at 30°C, washed with dH₂0 and resuspended in fresh YPD at a concentration of 1 x 10⁷ cells/ml. Cultures were subsequently arrested in late G1 phase with the addition of 5 μg/ml alpha-factor peptide (Sigma-Aldrich) for 2.5 hours. Cells were monitored for G1 phase arrest through microscope observation of the percentage of unbudded cells. Approximately 1.5 × 10⁷ cells were collected and treated with 0.1% sodium azide then kept on ice at 4°C. The remainder of the arrested culture was then centrifuged at 200 g and the pellet was washed twice with dH₂0. The cells were all BAR ⁺ and thus secrete the Bar1 protein, which degrades the alpha-factor pheromone used. Medium was collected and saved during centrifugation of the original logarithmic culture containing this protein and was subsequently used to resuspend the cells for the release from G1 phase. Additionally, Pronase E (Sigma-Aldrich), an enzyme that also hydrolyzes alpha-factor, was added at a concentration of 10 μg/ml to facilitate a synchronous release. Roughly the same number of cells collected in the G1-arrested samples were collected at the other time points following release and similarly treated with sodium azide and kept on ice until the completion of the time course.

Fluorescence activated cell sorting (FACS)

To assess cell synchrony, 1.5 × 10⁶ cells were removed from each time point sample. They were immediately centrifuged, re-suspended in 1 ml of ice-cold 70% ethanol and stored overnight at 4°C. Cells were then re-suspended in 500 μl of 50 mM Tris-HCl pH 8 containing 10 mg/ml of RNase A and incubated for 2 hours at 37°C. This was followed by centrifugation and re-suspension in 500 μl 50 mM Tris-HCl 7.5 with 2 mg/ml Proteinase K. Incubation at 50°C for an hour was performed prior to final resuspension in 100 μl FACS buffer (200 mM Tris-HCl 7.5, 200 mM NaCl and 78 mM MgCl₂). Cells were stained with SYTOX Green dye (5 μM; Molecular Probes) for at least one hour and then analyzed using a Becton-Dickinson FACScan.

Chromatin fractionation

Chromatin fractionation was performed as described by Semple et al. [8] with some modifications. Approximately 1 × 10⁷ cells collected from each time point were incubated in 7.5 ml pre- spheroplasting buffer (100 mM EDTA-KOH pH 8, 10 mM DTT), after washing once with dH₂O. They were then incubated at 30°C for 10 min with gentle shaking. Cells were centrifuged and re-suspended in 7.5 ml spheroplasting buffer (1 X YPD, 1.1 M sorbitol) containing 0.5 mg/ml Zymolyase 20 T (Seikagaku Corp., Japan) and 0.1 mg/ml Oxalyticase (Sigma), followed by shaking at 30°C for 30-45 min with gentle mixing. Cells were then washed once with 20 ml spheroplasting buffer containing 0.5 mM PMSF followed by resuspension in 1 ml ice-cold wash buffer (5 mM Tris–HCl pH 7.4, 20 mM KCl, 2 mM EDTA-KOH pH 7.4, 1 M sorbitol, 1% thiodiglycol, 125 mM spermidine, 50 mM spermine). Wash, Breakage and Lysis buffers all contained 1 tablet/10 ml of EDTA-free protease inhibitors (Roche) and were supplemented with 0.5 mM PMSF. Cells were centrifuged at 400 g for 1 min at 4°C, washed twice with 1 ml of Wash buffer and then re-suspended in 800 μl of Breakage buffer (5 mM Tris–HCl pH 7.4, 20 mM KCl, 2 mM EDTA-KOH pH 7.4, 0.4 M sorbitol, 1% thiodiglycol, 125 mM spermidine, 50 mM spermine). To these cells, 1 ml of Lysis buffer (Breakage buffer supplemented with 1% Triton X-100) was added and after repeated inversion (until the solution turned clear), cells were pelleted at 16,000 g for 10 min. This separated the proteins bound to the chromatin (residing in the pellet, referred to as PEL) from those solubilized in the non-chromatin fraction (supernatant or SUP). After removal of the supernatant, an additional 1 min spin at 16,000 g was performed to isolate any residual supernatant and the pellet was re-suspended in 100 μl Breakage buffer. MgCl₂ (5 mM) and DNase I (2 μg/ml) were added to the PEL fractions to solubilize the chromatin and associated proteins. After 10 min, the reaction was quenched with the addition of 2 μl 0.5 M EDTA. 10 μl of each SUP and PEL sample were collected for protein quantification. To the remaining of each sample, ½ volume of loading buffer (270 mM DTT, 9.9% SDS, 26% Glycerol and 10% Bromophenol blue) was added. Samples were subsequently boiled and equal concentrations of total protein were loaded on 7.5% SDS-PAGE gels. For an equal volume, PEL samples were 20-fold more concentrated than SUP samples due to the fact that approximately 5% of proteins are chromatin bound.

Western blotting and densitometry

Chromatin fractionation samples were assayed for protein concentration using the Bradford assay (Biorad). Equal amounts of protein were loaded into each lane of a 10% SDS-polyacrylamide gel for each set of time course samples. Following transfer, nitrocellulose membranes were stained with Ponceau S dye (Sigma). Membranes were destained using 1 X TEN buffer for 10 min. Mouse anti-Myc (Sigma, 1:5000) and Alexa Fluor 488 goat anti-mouse IgG (Invitrogen, 1:3000) antibodies were used to detect Myc-tagged Cdc45 and Cdc6. Anti-Mcm2 antibody (yN-19 goat polyclonal, Santa Cruz, 1:500) along with Alexa Fluor 488 donkey anti-goat IgG (Invitrogen, 1:3000) antibodies were used to detect Mcm2. Blots were also probed with TAT1 antibody (gift from the Gull lab) for visualizing α-tubulin, which should be exclusive to SUP samples as it does not associate with chromatin. Blots were incubated in primary and secondary antibodies for 2 h each, proceeding 4 h of blocking in 5% skim milk. Between blocking and each antibody treatment, blots were washed for 2 × 10 min with 1 x TEN + 0.05% Tween-20. A Typhoon 9410 scanner (GE healthcare) was used to analyze the blots. Densitometry readings were performed using Image Quant TL software (BD) and were normalized to total protein concentration as judged by the Ponceau S stain and/or tubulin band intensity. The same protocol for Western blotting was used in perturbation experiments involving HA-tagged strains. In this case, mouse anti-HA (Sigma, 1:5000) and Alexa Fluor 488 goat anti-mouse IgG (Invitrogen, 1:3000) antibodies were used to determine levels of Cdc6, Cdt1 and Dbf4.

Determination of in vivo cellular concentrations for each factor was performed for each time point. Normalization of densitometry readings was carried out by averaging the means of all SUP values from each set of experiments for a given protein. This average was divided by the means of the SUPs for each trial to give a scaling factor (S1) for each trial. Each SUP value for a given trial was then multiplied by its scaling factor (S1, S2 or S3). The same procedure was applied to all PEL sample values. In general, when a protein is in the chromatin fraction it is DNA-bound. To correct for non-specific DNA-binding, we determined, for each protein, a background level of non-specific binding corresponding to the observed abundance from a time-point at which the factor is known to be absent from origins. To obtain densitometric values, the program Imgage Quant™ was used to analyze blots scanned by a Typhoon™ 9400 imager (both GE Healthcare). To convert the densitometry measures to molecules/cell concentrations, we determined a scaling factor for each protein from the molecule counts reported in [67]. For each experiment, we averaged the total densitometry measure (SUP plus PEL) across the time-points to arrive at an averaged asynchronous densitometry reading (weighted according to the time contribution of each sample to a 90 minute cycle), which was then compared to the database to arrive at a scaling factor. We could not follow this procedure for Cdc6 or Dbf4, since they are not included in the Ghaemmmaghami et al. [67] database. A scaling factor was determined for these proteins by determining the relative abundance between Cdc45, Cdc6 and Dbf4 using similarly tagged asynchronous cultures. Whole cell extract levels of Cdc45-Myc and Cdc6-Myc were analysed by Western blotting and densitometry to yield an approximately 3:1 ratio of Cdc45 to Cdc6 copies per cell. Similarly, Cdc6-HA and Dbf4-HA were compared. This analysis allowed us to arrive at an asynchronous value of 576 molecules/cell for Cdc6 (Cdc6_Total = 576), i.e. one third of the value reported for the concentration of Cdc45, i.e. 1730 copies/cell according to [67] (Cdc45_Total = 1730). Similarly, the concentration of Dbf4 in an asynchronous population was 270 molecules/cell (Dbf4_Total = 270). The copy/cell number used for Mcm2 was 40,000 [58] while that for Cdt1 was 2190 [67] (Cdt1_Total = 2190). As an example of the agreement of our simulated values with literature-observed origin stoichiometry, Mcm2 (representing the MCM complex) was present at levels that were consistent with having two MCM complexes bound to each origin. This reflects the head to head placement of the heterododecamer at the origin. Once firing occurs, two forks are produced, illustrated in our model as one Mcm2-7·Cdc45 species molecule being generated (each FORK in the model represents a pair of these complexes). The level of chromatin-associated MCM protein that we obtained corresponds to roughly 300 origins being bound in this manner, which is the range of the number of origins that are reported to potentially fire per cell per cell cycle according to various global origin characterization studies [65, 90–92].

The values used to assign copies/cell numbers were taken from the best literature source available. The number for Mcm proteins in the GFP-tag database [67] differ by orders of magnitude from the number that has been reported by [58]. The latter is widely accepted, and so this discrepancy cannot be ignored. While similar inconsistencies might exist for other proteins in our model, where possible, we have carefully evaluated whether the GFP-tag database number gives a plausible cellular abundance through comparison with other reported values. For example, levels for Cdc45 in [67] are consistent with the value reported by [56], while no value has been reported for Cdt1 except for that in the GFP-tag database. In the event that future studies provide more accurate values for the various cellular abundances, the model can be easily revised to accommodate these changes by scaling the parameters associated with the relevant species. Indeed, it is common practice for models to be presented in terms of ‘arbitrary’ concentration units. Our use of a molecule-per-cell concentration scale allows us to make absolute predictions about species time-courses. Any future adjustments to the species levels will scale the concentrations and parameter estimates, but will have no direct bearing on the overall dynamic behaviour of the network.

Perturbation experiments

DY-139 (GAL1-CDC6), DY-140 (GAL1-CDT1) and DY-255 (GAL1-DBF4) strains were grown to 1 × 10⁷ cells/ml in 2% galactose/1% raffinose medium at 30°C, centrifuged at 6000 g for 5 min and washed with dH₂0. Cells were then resuspended in 2% glucose medium, maintaining the same cell concentration. Culture aliquots were removed for FACS analysis and preparation of whole cell extracts (as described in [93]), both before, and at various intervals after the switch to glucose. As references for normal endogenous protein levels, whole cell extracts were also prepared from strains expressing DBF4 (DY-256), CDC6 (DY-142) and CDT1 (DY-143) from their endogenous promoters. All strains produced fusion proteins with a 3HA epitope tag to facilitate visualization via Western blotting.

Time-varying model inputs

We collected protein data from yeast strains which were observed to have a generation time of ~90 min in log phase. In order to fit our data to the timescale of the Chen et al. [45] model, which has a period of 101.2 minutes, we scaled our experimental time-course to reflect this change in timing. This value, as described in their model, is chosen to reflect the longer cycle of daughter cells (which are smaller than mother cells in asymmetric cell division). Time-point samples were collected from an alpha factor arrest in late G1 phase, corresponding to experimental time-point T = 0. Additional samples were collected at 5, 10, 15, 30, 45, 60 and 75 min after synchronous release from the alpha factor block. We observed that our experimental timepoint T = 0 corresponds to 19 min after the beginning of G1 phase in the 101.2 min model. We arrived at this number by comparing the point at which cells entered S phase in vivo (~15 min after alpha factor release as determined by FACS analysis) and the corresponding start of S phase in the model.

The difference between the original Cdc20 profile and our modified version is shown in Additional file 8.

Implementation of a combined model

In order to implement a combined model, we closed the loop between our model and the cell cycle model of Chen et al. [45] by eliminating the lumped species ORI and replacing it with an indication of the progress of replication from our FORK species. In this combined model we modified the formula for k_mad2 shown above by replacing the ratio (ORI/15) with $(\int FORK/500)$, where the integration begins at the start of the cell cycle. Since each model incorporates a description of Cdc6 dynamics, we merged them by including the dynamics of origin binding from our model with the Cdc6 dynamics of the Chen model. This resulted in good accordance between the behaviour of the replication initiation network in our model in the combined model. There are two marginal differences, neither of which affects the overall system behaviour (see Additional file 8). Firstly, due to the fact that our internal model uses a Clb5 profile generated by scaling from arbitrary units, it is nearly identical to the Chen model Clb5 profile. The Chen Clb5 profile extends farther past the 101.2 minute mark than ours resulting in a more rapid dephosphorylation of RC7 in the internal model. Thus, RC7 persists for ~5 mins longer in the combined model. Secondly, because the merged Cdc6 decreases in concentration earlier than in our model in isolation, RC1 levels stay high until fork firing occurs. Despite these differences, the essential dynamics of the system are preserved: replication fork firing follows the same pattern, with the RC7→RC1 delay having no effect on timing of firing. The behaviours of the Chen model species are not perceptibly altered by the removal of ORI and the combination of the two models (Additional file 6).

The effects of combining our models are shown to be minimal in the wildtype case. Comparing eight cell cycle mutants used to fit the Chen model in the combined model, we see the identical phenotypes (see Additional file 7 and 9). These results demonstrate that our replication initiation module is a functional replacement for the corresponding black box in Chen’s whole cell cycle model. A summary of simulated mutant phenotypes [7, 55, 60, 73, 74, 94–97] is given in Additional file 9.