Stochastic hybrid modeling of DNA replication across a complete genome

A Stochastic Hybrid Model of DNA Replication.

Stochastic hybrid systems are used in areas such as mathematical finance, transportation, and telecommunications to model systems that involve the interaction of discrete phenomena, continuous phenomena, and uncertainty (28). DNA replication is a stochastic hybrid process, where discrete transitions between origin states are coupled with the continuous movement of replication forks along the genome, while initiation events are characterized by uncertainty in time and space. We have built a stochastic hybrid model that reproduces the DNA replication process genome-wide. The basic features of the model are presented here; a full description including formal definitions is given in supporting information (SI) Text, see also (28) and (29) for the mathematical foundations on which the model is based.

The discrete dynamics of the model capture instantaneous changes in the state of the system, such as origin firing or fork conversion. The state of each origin i is represented by a variable, S_i, that takes one of six values (PreR, PassR, RB, RL, RR, PostR) (Fig. 1). All origins start at the Pre-Replicative state (PreR) and end either in the passively replicated (PassR) or the Post-Replicating (PostR) states. The remaining states discriminate origins from which active forks emanate on both directions (RB), only to the left (LR) or only to the right (RL). Transitions between states are governed by guards (G), logical statements that depend on the continuous and stochastic dynamics of the system and determine when the transition will take place. The continuous dynamics of the model capture the evolution of replication forks along the genome and are represented by two differential equations for each active origin, denoting the number of bases replicated by forks moving away from origin i toward the left (L_i) or right (R_i). Notice that there is a tight coupling between the discrete and continuous states of the process, a defining feature of hybrid systems. Origin firing is characterized by a degree of uncertainty in both the location and time of firing of different origins, captured by the stochastic dynamics of the model. For example, the time T_i at which each origin i will fire in any given cell in a population can be randomly extracted based on its experimentally defined firing probability. Note that each origin is characterized by an intrinsic propensity to fire (referred to hereafter as firing propensity), which denotes the instantaneous probability of firing in the absence of passive replication. The firing propensity differs from the observed firing efficiency (referred to hereafter as observed efficiency), which also depends on the probability of passive replication for this particular origin.

Model Flexibility.

The model developed is generally applicable and can, in principle, capture the DNA replication process of any organism for which knowledge of the process of origin selection and activation is available. The model requires as input:

1. The length of the genome in base pairs (L).

2. The number (N) of putative origins of replication and the location X_i of each putative origin i along the genome. The number and location of putative origins can be either deterministic based on experimental data, or extracted randomly based on information such as the average number or average spacing of origins.

3. The firing propensity λ_i of each origin i. This is related to the probability of the origin firing in a unit of time in the absence of passive replication. The model can accommodate different modes of origin firing in time and space (such as stochasticity in space versus a spacing mechanism, stochasticity in time versus a timing mechanism, or constant firing propensity during S-phase versus increasing or decreasing firing propensity as a function of time). We will provide an example of this below.

4. The speed v(x) with which replication forks progress when they go over position x of the genome. Speed can be deterministically defined, or depend on system parameters, such as interorigin spacing.

Model Instantiation.

We have used experimental data from the fission yeast Schizosaccharomyces pombe for model instantiation. In this organism, there is experimental support for a stochastic mode of origin activation provided at the single cell level (27). High-throughput analyses have permitted the location of the majority of putative origins to be specified along the genome (23–25). Additionally, the firing propensity of each origin along the genome has been determined based on the experimentally defined observed efficiencies in the presence of the drug hydroxyurea, which inhibits fork progression, thereby eliminating passive replication (24). In fission yeast, only a small number of origins genome-wide appear intrinsically late and do not fire during a long hydroxyurea block (24, 27), while for the vast majority of origins, the timing of firing appears correlated to their firing propensities (24), suggesting that stochastic determination of firing times is the prevailing mechanism in the fission yeast (18). Fission yeast offers therefore a good opportunity to capture full genome replication based on experimental data coupled to an understanding of the underlining biology.

The number of origins (N = 893), their exact locations (X₁, X₂, …, X_N) along the three chromosomes of S. pombe (total length L = 12,039,987 bases), and the firing propensity of each origin (FP_i, for i=1, 2, …, N, defined as the fraction of cells in which each origin was observed to fire in hydroxyurea) were used as input (24), see Table S1. A constant speed of v(x) = 3000 bases/minute was used genome-wide, based on experimental estimates (24). The firing time, T_i (in minutes), of each origin i was extracted according to an exponential distribution whose rate λ_i depends on FP_i (SI Text). Initially, the firing times of different origins in the absence of passive replication were assumed to be statistically independent. This is a valid assumption for fission yeast, based on available experimental data (24). In addition, the firing propensity of each origin as a function of time (λ_i) was assumed to be constant during S-phase. This means that the total remaining firing propensity of the system decreases as S-phase progresses. An alternative model that relaxes both assumptions will be discussed below (Propensity Redistribution Model).

Simulation Results.

The model was coded in the Matlab environment. Monte-Carlo simulations were used to capture the process of DNA replication at a single cell, extract useful statistics on model predictions at a population level, and validate them against experimental data. Examples of simulation runs are given in Fig. S1 and show the stochastic nature of the process. Diagnostic tests were carried out which showed that observed fluctuations in continuous and discrete system parameters during the evolution of DNA replication behave as predicted and each part of the genome is copied once and only once (see SI Text and Fig. S2).

To validate the model using an independent source of experimental data, the output of simulations was compared to data derived by the analysis of replication intermediated at the single molecule level by DNA combing (27). In this method, replication intermediates are labeled in vivo, genomic DNA is extracted, DNA fibers stretched and visualized. To reproduce in silico the experimental conditions of (27), random genomic pieces with a size distribution matching the one of the fibers analyzed in (27) were extracted from 2,000 simulations under three circumstances: assuming that the first 160 origins (24) scheduled to fire in each simulation would arrest after 5 kb [mimicking the hydroxyurea block of (27)], when 30% of the genome had been replicated in each simulation (early S-phase) and when 70% of the genome had been replicated in each simulation (late S-phase). In Fig. 2, the distribution of distances between replication bubbles (inter-bubble distances) is shown for each case. The exponential curves generated by (27) fitting their experimental data have been overlaid for comparison. In all cases, simulated results closely matched the experimentally derived curves (R² values of 0.9937, 0.9794, and 0.9749, respectively).

Model output was used to analyze genome-wide replication kinetics in silico. In Fig. 3A, the fraction of unreplicated DNA as a function of time is shown, with each curve representing one simulation (single cell). While in all simulations the process of conversion from fully unreplicated (1) to fully replicated (0) follows similar kinetics, differences in individual cells due to stochastic phenomena are evident. Strikingly, while the bulk of DNA replication takes place quickly during the first half of the process, a very long tail is evident in the second half, showing that completion of DNA replication of a remaining small fraction of the genome takes a long time. This is also evident in the rate of new DNA synthesis (Fig. 3B), which shows the expected dumb-bell shape observed experimentally, but continues with a long tail. As a consequence, the mean completion time of the process (Fig. 3C) is unexpectedly long [67 min, as opposed to experimental estimates of ≈20 min (30)]. This is not due to less origins firing in each simulation than expected (Fig. 3D, average 159 origins firing in each simulation, very close to the experimental estimate of ≈160). In each simulation, at least one region of the genome exhibits a distance between adjacent active origins that is very long (Fig. 3E, mean maximum interorigin distance of 258 kb, replication of which would last 43 min with a fork speed of 3,000 bases/minute, assuming both origins fire at the same time). Different locations along each chromosome delay the process in different simulations (Fig. 3F, showing regions that terminate replication after 99% of the genome has replicated) and even the most problematic regions do not delay the process in >12% of the simulations. This is explained by the random distribution of active origins along the genome (random gap problem). The model predicts that there is a small percentage (<5%) of DNA that remains unreplicated for a long time, that the rate of DNA synthesis at late time points is low and that there is an infrequent appearance within the population of each problematic region. Such a delay in S-phase completion is difficult to detect experimentally and could well have gone unnoticed.

Sensitivity Analysis.

Given the unexpectedly long S-phase completion times revealed by full genome simulations, we proceeded to vary the input parameters to test their effect on system behavior. Fig. 4 shows how the predicted S-phase duration and the number of firing origins depends on three of the model parameters: the intrinsic firing propensity of each origin; the fork speed v(x); and the number of origins N.

To systematically increase the firing propensity of all origins we decreased the parameter T_f (a system parameter inversely correlated to the firing propensity of each origin, see SI Text) from 20 min down to 1 min (Fig. 4 A–C). As firing propensities increase, the predicted average S-phase duration decreases, from 67 to 21 min. As expected, however, the predicted average number of firing origins also increases from 18% to 57%. While small values of T_f (high firing propensities) lead to reasonable S-phase completion times, the high number of firing origins predicted for these values is not supported by experimental data (22, 24, 27, 31).

As fork speed v(x) increases from 1 kb/min to ≈6 kb/min, S-phase completion times are quickly reduced (from an average of 126 min to an average of 44 min, Fig. 4D). However, the fraction of firing origins is similarly reduced (Fig. 4E); moreover, further increases in speed have smaller effects on the overall completion time. For S-phase completion time to reach the expected 20 min, speeds of >18 kb/min would be required, which are unrealistically high.

Finally, the number of potential origins, N, was increased by introducing more low efficiency origins at specific locations along the genome (Fig. 4 G–I). The experimental method used for origin identification (24) could not reliably detect origins with <10% efficiency. New origins with efficiencies below this experimental threshold were introduced into biologically relevant chromosomal locations, predicted by bioinformatics analysis to fulfill characteristics of known origin regions (SI Text). Criteria depending on AT content and size were used to gradually introduce additional origins with observed efficiencies of FP_i = 8%. The addition of new origins causes a decrease in the predicted average S-phase duration, however even after adding all potential origins predicted by our analysis (≈50% of all intergenes in the fission yeast genome), S-phase completion times still remained long (46 min).

Firing Propensity Redistribution Model.

Our analysis leads to two alternatives: Either S-phase duration in the presence of stochastic firing is longer than experimentally defined, or a further mechanism exists that succeeds in limiting S-phase duration circumventing the random gap problem. One proposed solution to the random gap problem is that the efficiency of unreplicated origins may increase during the course of S-phase (14, 18). Stochasticity could be brought about by a limiting factor that allows firing of a fraction of putative origins. The redistribution of this factor after an origin fires or gets passively replicated would imply a gradual increase in firing propensity of origins still in the prereplicative state as S-phase progresses (18). This biological scenario can be captured in our model by keeping the total system firing propensity (given by the sum of the firing propensities of all origins) constant throughout S-phase. Each origin that fires or becomes passively replicated releases its firing propensity λ_i, which is redistributed to the remaining PreR origins proportionally to their initial firing propensities. We call the resulting variant of the model the “Firing Propensity Redistribution Model” and we describe it in detail in the SI Text.

Fig. 5A shows the evolution of unreplicated DNA over time predicted by the firing propensity redistribution model. While the kinetics of bulk DNA replication are similar to the basic model, the curve has lost its long tail. A similar change is observed in the curve depicting the rate of newly synthesized DNA over time (Fig. 5B). Consequently, the average time of completion of DNA replication has decreased from 67 min in the basic model to 33 min (Fig. 5C), concomitant with a decrease in the mean maximum distance observed between active origins in each simulation (from 258 kb to 135 kb, Fig. 5E). The number of firing origins has increased to a mean of 286 (Fig. 5D). The regions that replicate last in the process are shown in Fig. 5F. Two additional variants of the redistribution model (redistribution upon fork conversion proportional to initial firing propensity and redistribution upon firing proportional to current firing propensity) are discussed in the SI Text. The simulation results for the two variants (data not shown) were similar to the results of Fig. 5, with redistribution upon fork conversion leading to slightly longer S-phase duration (average replication completion time of 37 min).

We conclude that the redistribution during S-phase of a limiting factor acting at the level of initiation of DNA replication is one way by which the problems arising due to stochastic origin activation can be circumvented.

Simulating Full Genome DNA Replication.

A general and flexible stochastic hybrid model was developed to capture the DNA replication process during the S-phase of the cell cycle. The model in principle applies to any organism; here it was instantiated using experimental data for Schizosaccharomyces pombe (24). The resulting model was coded in simulation and tested against independent genome-wide data on fission yeast DNA replication (27), showing that it captures well our current understanding of the DNA replication process. The predictions of the model suggest that given stochastic firing, S-phase duration would be considerably longer than generally accepted for this organism [on average 67 min and reaching 113 min in some cells, instead of 20 min (30)], due to randomly generated large interorigin gaps.

The Effect of Parameter Values.

The model's flexibility allowed us to introduce modifications to test possible explanations for this difference. Input parameters (intrinsic firing propensities of all origins, fork speed and number of potential origins) were varied to investigate their effect on S-phase duration. While the parameters tested have a direct effect on S-phase duration, none of these changes alone seems capable of explaining in a biologically realistic way the discrepancy between experimentally measured S-phase duration and the predictions of the model. Local effects, such as an increase in fork speed at problematic regions, could conceivably speed up DNA replication completion. However, different regions create problems in different cells (see Fig. 3F) and local effects at a number of these problematic regions would be required for an effect at the population level. Combining changes in input parameters could potentially reveal a biologically plausible solution, although we have been unable to pinpoint a likely combination (data not shown). More complex patterns of origin firing, such as coordination of fork velocity with interorigin distance (17) or defined spacing of active origins (15, 16) could potentially decrease S-phase duration. Such conjectures can be tested by our model and we have made use of the model's versatility to capture one such possibility.

Firing Propensity Redistribution.

The model was modified to accommodate a specific biological hypothesis for the reason behind the probabilistic nature of the firing of origins: the presence of a limiting initiation factor, which is released following origin firing and binds again to origins still at the prereplicative state (18). Since the number of origins in the prereplicative state decreases as S-phase progresses, the probability of the limiting factor binding to a particular prereplicative origin (and therefore the firing probability of this origin) increases along S-phase (18). The firing propensity redistribution model was developed to code this hypothesis, assuming that the pool of available factor remains roughly constant during S-phase. Using the initial model parameters, the firing propensity redistribution model predicts a considerable reduction in S-phase duration, while also maintaining the average fraction of origins that fire during S-phase in a reasonable range. Strikingly, while bulk replication progresses with dynamics similar to the basic model, the long tail of replication has disappeared. The firing propensity redistribution model reproduces the experimental data of (27) as well as the base-case model does (data not shown), consistent with the bulk of replication progressing in a similar manner in the two models. Firing propensity redistribution therefore offers a possible solution to the random completion problem, without the need to postulate a spacing mechanism for positioning active origins at appropriate locations. This scenario is indeed biologically plausible. The limiting factor should act at the level of origin activation, and not origin licensing that is globally inhibited following S-phase onset. It could be released at firing (Fig. 5) or upon fork conversion (data not shown). Cyclin-dependent kinases or the Dbf4-dependent kinase (DDK), which are required for initiation, are possible candidates for the first option, while Cdc45, previously shown to be limiting for replication in Xenopus (32), is a candidate for the second option.

“G2” Replication?

An alternative explanation is that DNA replication takes longer than has been generally accepted. The predicted low percentage of unreplicated DNA (<5%, Fig. 3A), the very low rate of replication at late time points (Fig. 3B), and the fact that problematic regions (those that delay the overall process) differ in different cells (Fig. 3F), suggest that such a prolonged phase of DNA replication could have gone unnoticed. The inability of techniques used up to now to detect small amounts of unreplicated DNA scattered across different regions in different cells means that the length of S-phase may be longer than our current estimates.

An S-phase that extends into what we now define as G2-phase is consistent with known properties of the fission yeast cell cycle and the available experimental data (Fig. 6). G2 makes up the largest part of the fission yeast cell cycle, and it would be inconsequential if small parts of the genome continued replication for an extended time while the cell grows in size. The checkpoint that arrests entry into mitosis in the presence of unreplicated DNA would become essential if the G2/M size control was inoperative. This alternative explanation can explain a number of previous observations in fission yeast. In fission yeast, S+G2 is not shortened beneath a certain minimal period (0.5 of the cell cycle) even in large cells that have reached the critical size for mitosis (33, 34). Furthermore, a wee1 mutant in which the G2/M size control is inoperative and, as a consequence, advances into mitosis at a small size (33, 35) also has a S+G2 length of approximately 0.5 of the cell cycle. Significantly, wee1 is synthetically lethal with a rad3 mutant, defective in the G2/M checkpoint control (36). We postulate that in wee1 mutant cells, the extended S-phase occupies all of the previously designated “G2” period and mitosis is restrained by the Rad3-dependent G2/M checkpoint. When the checkpoint is inactivated, cells cannot restrain entry into mitosis even though S is not complete, and this leads to lethality. The minimal S+G2 (65–120 min depending on growth conditions), is similar to the extended S-phase suggested here (Fig. 3A). Additional support for the extended S-phase hypothesis comes from studies of New-End Take-Off (NETO), a growth transition that is dependent on the completion of DNA replication and takes place in early G2. This occurs at approximately 0.35 of the cycle (37) rather than 0.1 of the cycle, the time when DNA replication is believed to be completed (30). NETO shows heterogeneity among cells in a population (38), consistent with the results presented here.

Any organism with an active checkpoint inhibiting entry into mitosis until DNA replication has been completed might not require a further specific mechanism to inhibit the random gap problem. In this light, part or all of what we now experimentally define as G2-phase might be better described as the period required for the cell to complete replication of the last remaining parts of its genome. Proving or disproving the extended S-phase hypothesis will have to await the development of more sensitive methods, able to detect the low amounts of late replication suggested by the model described here.