The eukaryotic bellshaped temporal rate of DNA replication origin firing emanates from a balance between origin activation and passivation
Abstract
The timedependent rate $I\left(t\right)$ of origin firing per length of unreplicated DNA presents a universal bell shape in eukaryotes that has been interpreted as the result of a complex timeevolving interaction between origins and limiting firing factors. Here, we show that a normal diffusion of replication fork components towards localized potential replication origins (poris) can more simply account for the $I\left(t\right)$ universal bell shape, as a consequence of a competition between the origin firing time and the time needed to replicate DNA separating two neighboring poris. We predict the $I\left(t\right)$ maximal value to be the product of the replication fork speed with the squared pori density. We show that this relation is robustly observed in simulations and in experimental data for several eukaryotes. Our work underlines that forkcomponent recycling and potential origins localization are sufficient spatial ingredients to explain the universality of DNA replication kinetics.
https://doi.org/10.7554/eLife.35192.001eLife digest
Before a cell can divide, it must duplicate its DNA. In eukaryotes – organisms such as animals and fungi, which store their DNA in the cell’s nucleus – DNA replication starts at specific sites in the genome called replication origins. At each origin sits a protein complex that will activate when it randomly captures an activating protein that diffuses within the nucleus. Once a replication origin activates or “fires”, the complex then splits into two new complexes that move away from each other as they duplicate the DNA. If an active complex collides with an inactive one at another origin, the latter is inactivated – a phenomenon known as origin passivation. When two active complexes meet, they release the activating proteins, which diffuse away and eventually activate other origins in unreplicated DNA.
The number of origins that activate each minute divided by the length of unreplicated DNA is referred to as the “rate of origin firing”. In all eukaryotes, this rate – also known as I(t) – follows the same pattern. First, it increases until more than half of the DNA is duplicated. Then it decreases until everything is duplicated. This means that, if plotted out, the graph of origin firing rate would always be a bellshaped curve, even for organisms with genomes of different sizes that have different numbers of origins. The reason for this universal shape remained unclear.
Scientists had tried to create numerical simulations that model the rate of origin firing. However, for these simulations to reproduce the bellshape curve, a number of untested assumptions had to be made about how DNA replication takes place. In addition, these models ignored the fact that it takes time to replicate the DNA between origins.
To take this time into account, Arbona et al. instead decided to model the replication origins as discrete and distinct entities. This way of building the mathematical model succeeded in reproducing the universal bell curve shape without additional assumptions. With this simulation, the balance between origin activation and passivation is enough to achieve the observed pattern.
The new model also predicts that the maximum rate of origin firing is determined by the speed of DNA replication and the density of origins in the genome. Arbona et al. verified this prediction in yeast, fly, frog and human cells – organisms with different sized genomes that take between 20 minutes and 8 hours to replicate their DNA. Lastly, the prediction also held true in yeast treated with hydroxyurea, an anticancer drug that slows DNA replication.
A better understanding of DNA replication can help scientists to understand how this process is perturbed in cancers and how drugs that target DNA replication can treat these diseases. Future work will explore how the 3D organization of the genome affects the diffusion of activating proteins within the cell nucleus.
https://doi.org/10.7554/eLife.35192.002Introduction
Eukaryotic DNA replication is a stochastic process (Hyrien et al., 2013; Hawkins et al., 2013; Hyrien, 2016b). Prior to entering the S(ynthesis)phase of the cell cycle, a number of DNA loci called potential origins (poris) are licensed for DNA replication initiation (Machida et al., 2005; Hyrien et al., 2013; Hawkins et al., 2013). During Sphase, in response to the presence of origin firing factors, pairs of replication forks performing bidirectional DNA synthesis will start from a subset of the poris, the active replication origins for that cell cycle (Machida et al., 2005; Hyrien et al., 2013; Hawkins et al., 2013). Note that the inactivation of poris by the passing of a replication fork called origin passivation, forbids origin firing in already replicated regions (de Moura et al., 2010; Hyrien and Goldar, 2010; Yang et al., 2010). The timedependent rate of origin firing per length of unreplicated DNA, $I\left(t\right)$, is a fundamental parameter of DNA replication kinetics. $I\left(t\right)$ curves present a universal bell shape in eukaryotes (Goldar et al., 2009), increasing toward a maximum after midSphase and decreasing to zero at the end of Sphase. An increasing $I\left(t\right)$ results in a tight dispersion of replication ending times, which provides a solution to the random completion problem (Hyrien et al., 2003; Bechhoefer and Marshall, 2007; Yang and Bechhoefer, 2008).
Models of replication in Xenopus embryo (Goldar et al., 2008; Gauthier and Bechhoefer, 2009) proposed that the initial $I\left(t\right)$ increase reflects the progressive import during Sphase of a limiting origin firing factor and its recycling after release upon forks merge. The $I\left(t\right)$ increase was also reproduced in a simulation of human genome replication timing that used a constant number of firing factors having an increasing reactivity through Sphase (Gindin et al., 2014). In these three models, an additional mechanism was required to explain the final $I\left(t\right)$ decrease by either a subdiffusive motion of the firing factor (Gauthier and Bechhoefer, 2009), a dependency of firing factors’ affinity for poris on replication fork density (Goldar et al., 2008), or an inhomogeneous firing probability profile (Gindin et al., 2014). Here, we show that when taking into account that poris are distributed at a finite number of localized sites then it is possible to reproduce the universal bell shape of the $I\left(t\right)$ curves without any additional hypotheses than recycling of fork components. $I\left(t\right)$ increases following an increase of fork mergers, each merger releasing a firing factor that was trapped on DNA. Then $I\left(t\right)$ decreases due to a competition between the time ${t}_{c}$ to fire an origin and the time ${t}_{r}$ to replicate DNA separating two neighboring pori. We will show that when ${t}_{c}$ becomes smaller than ${t}_{r}$, pori density over unreplicated DNA decreases, and so does $I\left(t\right)$. Modeling random localization of active origins in Xenopus embryo by assuming that every site is a (weak) pori, previous work implicitly assumed ${t}_{r}$ to be close to zero (Goldar et al., 2008; Gauthier and Bechhoefer, 2009) forbidding the observation of a decreasing $I\left(t\right)$. Licensing of a limited number of sites as pori thus appears to be a critical property contributing to the observed canceling of $I\left(t\right)$ at the end of Sphase in all studied eukaryotes.
Results
Emergence of a bellshaped $I\left(t\right)$
In our modeling of replication kinetics, a bimolecular reaction between a diffusing firing factor and a pori results in an origin firing event; then each half of the diffusing element is trapped and travels with a replication fork until two converging forks merge (termination, Figure 1a). A molecular mechanism explaining the synchronous recruitment of firing factors to both replication forks was recently proposed (Araki, 2016), supporting the bimolecular scenario for pori activation. Under the assumption of a wellmixed system, for every time step $dt$, we consider each interaction between the ${N}_{FD}\left(t\right)$ free diffusing firing factors and the ${N}_{\text{p}\text{ori}}\left(t\right)$ poris as potentially leading to a firing with a probability ${k}_{on}dt$. The resulting simulated firing rate per length of unreplicated DNA is then:
where ${N}_{fired}(t,t+dt)$ is the number of poris fired between times $t$ and $t+dt$, and ${L}_{unrepDNA}\left(t\right)$ is the length of unreplicated DNA a time $t$. Then we propagate the forks along the chromosome with a constant speed $v$, and if two forks meet, the two half firing complexes are released and rapidly reform an active firing factor. Finally, we simulate the chromosomes as 1D chains where prior to entering Sphase, the poris are precisely localized. For Xenopus embryo, the pori positions are randomly sampled, so that each simulated Sphase corresponds to a different positioning of the poris. We compare results obtained with periodic or uniform pori distributions (Materials and methods). For S. cerevisiae, the pori positions, identical for each simulation, are taken from the OriDB database (Siow et al., 2012). As previously simulated in human (Löb et al., 2016), we model the entry in Sphase using an exponentially relaxed loading of the firing factors with a time scale shorter than the Sphase duration ${T}_{phase}$ (3 min for Xenopus embryo, where ${T}_{phase}\sim 30$ min, and 10 min for S. cerevisiae, where ${T}_{phase}\sim 60$ mins). After the short loading time, the total number of firing factors ${N}_{D}^{T}$ is constant. As shown in Figure 1b (see also Figure 2), the universal bell shape of the $I\left(t\right)$ curves (Goldar et al., 2009) spontaneously emerges from our model when going from weak to strong interaction, and decreasing the number of firing factors below the number of poris. The details of the firing factor loading dynamics do not affect the emergence of a bell shaped $I\left(t\right)$, even though it can modulate its precise shape, especially early in Sphase.

Figure 2—source data 1
 https://doi.org/10.7554/eLife.35192.006

Figure 2—source data 2
 https://doi.org/10.7554/eLife.35192.007

Figure 2—source data 3
 https://doi.org/10.7554/eLife.35192.008
In a simple bimolecular context, the rate of origin firing is $i\left(t\right)={k}_{on}{N}_{\text{p}\text{ori}}\left(t\right){N}_{FD}\left(t\right)$. The firing rate by element of unreplicated DNA is then given by
where ${\rho}_{\text{p}\text{ori}}\left(t\right)={N}_{\text{p}\text{ori}}\left(t\right)/{L}_{unrepDNA}\left(t\right)$. In the case of a strong interaction and a limited number of firing factors, all the diffusing factors react rapidly after loading and ${N}_{FD}\left(t\right)$ is small (Figure 1 (c), dashed curves). Then follows a stationary phase where as long as the number of poris is high (Figure 1 (c), solid curves), once a diffusing factor is released by the encounter of two forks, it reacts rapidly, and so ${N}_{FD}\left(t\right)$ stays small. Then, when the rate of fork mergers increases due to the fact that there are as many active forks but a smaller length of unreplicated DNA, the number of free firing factors increases up to ${N}_{D}^{T}$ at the end of Sphase. As a consequence, the contribution of ${N}_{FD}\left(t\right)$ to $I\left(t\right)$ in Equation (2) can only account for a monotonous increase during the S phase. For $I\left(t\right)$ to reach a maximum ${I}_{max}$ before the end of Sphase, we thus need that ${\rho}_{\text{p}\text{ori}}\left(t\right)$ decreases in the late Sphase. This happens if the time to fire a pori is shorter than the time to replicate a typical distance between two neighboring poris. The characteristic time to fire a pori is ${t}_{c}=1/{k}_{on}{N}_{FD}\left(t\right)$. The mean time for a fork to replicate DNA between two neighboring poris is ${t}_{r}=d\left(t\right)/v$, where $d\left(t\right)$ is the mean distance between unreplicated poris at time $t$. So the density of origins is constant as long as:
or
Thus, at the beginning of the Sphase, ${N}_{FD}\left(t\right)$ is small, ${\rho}_{\text{p}\text{ori}}\left(t\right)$ is constant (Figure 1 (c), solid curves) and so ${I}_{S}\left(t\right)$ stays small. When ${N}_{FD}\left(t\right)$ starts increasing, as long as Equation (4) stays valid, ${I}_{S}\left(t\right)$ keeps increasing. When ${N}_{FD}\left(t\right)$ becomes too large and exceeds ${N}_{FD}^{\ast}$, then Equation (4) is violated and the number of poris decreases at a higher rate than the length of unreplicated DNA, and ${\rho}_{\text{p}\text{ori}}\left(t\right)$ decreases and goes to zero (Figure 1 (c), red solid curve). As ${N}_{FD}\left(t\right)$ tends to ${N}_{D}^{T}$, ${I}_{S}\left(t\right)$ goes to zero, and its global behavior is a bell shape (Figure 1 (b), red). Let us note that if we decrease the interaction strength (${k}_{on}$), then the critical ${N}_{FD}^{\ast}$ will increase beyond ${N}_{D}^{T}$ (Figure 1 (c), dashed blue and green curves). ${I}_{S}\left(t\right)$ then monotonously increase to reach a plateau (Figure 1 (b), green), or if we decrease further ${k}_{on}$, ${I}_{S}\left(t\right)$ present a very slow increasing behavior during the Sphase (Figure 1 (b), blue). Now if we come back to strong interactions and increase the number of firing factors, almost all the poris are fired immediately and ${I}_{S}\left(t\right)$ drops to zero after firing the last pori.
Another way to look at the density of poris is to compute the ratio of the number of passivated origins by the number of activated origins (Figure 1 (d)). After the initial loading of firing factors, this ratio is higher than one. For weak and moderate interactions (Figure 1 (d), blue and green solid curves, respectively), this ratio stays bigger than one during all the Sphase, where ${I}_{S}\left(t\right)$ was shown to be monotonously increasing (Figure 1 (b)). For a strong interaction (Figure 1 (b), red solid curve), this ratio reaches a maximum and then decreases below one, at a time corresponding to the maximum observed in ${I}_{S}\left(t\right)$ (Figure 1 (d), red solid curve). Hence, the maximum of $I\left(t\right)$ corresponds to a switch of the balance between origin passivation and activation, the latter becoming predominant in late Sphase. We have seen that up to this maximum ${\rho}_{\text{p}\text{ori}}\left(t\right)\approx cte\approx {\rho}_{0}$, so ${I}_{S}\left(t\right)\approx {k}_{on}{\rho}_{0}{N}_{F}\left(t\right)$. When ${N}_{FD}\left(t\right)$ reaches ${N}_{FD}^{\ast}$, then ${I}_{S}\left(t\right)$ reaches its maximum value:
where we have used the approximation $d\left(t\right)\approx d\left(0\right)=1/{\rho}_{0}$ (which is exact for periodically distributed poris). ${I}_{max}$ can thus be predicted from two measurable parameters, providing a direct test of the model.
Comparison with different eukaryotes
Xenopus embryo
Given the huge size of Xenopus embryo chromosomes, to make the simulations more easily tractable, we rescaled the size $L$ of the chromosomes, ${k}_{on}$ and ${N}_{D}^{T}$ to keep the duration of Sphase ${T}_{phase}\approx L/2v{N}_{D}^{T}$ and $I\left(t\right)$ (Equation (2)) unchanged ($L\to \alpha L$, ${N}_{D}^{T}\to \alpha {N}_{D}^{T}$, ${k}_{on}\to {k}_{on}/\alpha $). In Figure 2 (a) are reported the results of our simulations for a chromosome length $L=3000$ kb. We see that a good agreement is obtained with experimental data (Goldar et al., 2009) when using either a uniform distribution of poris with a density ${\rho}_{0}=0.70$ kb${}^{1}$ and a number of firing factors ${N}_{D}^{T}=187$, or a periodic distribution with ${\rho}_{0}=0.28$ kb${}^{1}$ and ${N}_{D}^{T}=165$. A higher density of poris was needed for uniformly distributed poris where $d\left(t\right)$ (slightly) increases with time, than for periodically distributed poris where $d\left(t\right)$ fluctuates around a constant value $1/{\rho}_{0}$. The uniform distribution, which is the most natural to simulate Xenopus embryo replication, gives a density of activated origins of 0.17 kb${}^{1}$ in good agreement with DNA combing data analysis (Herrick et al., 2002) but twice lower than estimated from real time replication imaging of surfaceimmobilized DNA in a soluble Xenopus egg extract system (Loveland et al., 2012). Note that in the latter work, origin licensing was performed in condition of incomplete chromatinization and replication initiation was obtained using a nucleoplasmic extract system with strong initiation activity, which may explain the higher density of activated origins observed in this work.
S. cerevisiae
To test the robustness of our minimal model with respect to the distribution of poris, we simulated the replication in S. cerevisiae, whose poris are known to be well positioned as reported in OriDB (Siow et al., 2012). 829 poris were experimentally identified and classified into three categories: Confirmed origins (410), Likely origins (216), and Dubious origins (203). When comparing the results obtained with our model to the experimental $I\left(t\right)$ data (Goldar et al., 2009) (Figure 2 (b)), we see that to obtain a good agreement we need to consider not only the Confirmed origins but also the Likely and the Dubious origins. This shows that in the context of our model, the number of poris required to reproduce the experimental $I\left(t\right)$ curve in S. cerevisiae exceeds the number of Confirmed and Likely origins. Apart from the unexpected activity of Dubious origins, the requirement for a larger number of origins can be met by some level of random initiation (Czajkowsky et al., 2008) or initiation events away from mapped origins due to helicase mobility (Gros et al., 2015; Hyrien, 2016a). If fork progression can push helicases along chromosomes instead of simply passivating them, there will be initiation events just ahead of progressing forks. Such events are not detectable by the replication profiling experiments used to determine $I\left(t\right)$ in Figure 2(b) and thus not accounted for by ${I}_{max}$. Given the uncertainty in replication fork velocity (a higher fork speed would require only Confirmed and Likely origins) and the possible experimental contribution of the poris in the rDNA part of chromosome 12 (not taken into account in our modeling), this conclusion needs to be confirmed in future experiments. It is to be noted that even if 829 poris are needed, on average only 352 origins have fired by the end of Sphase. For S. cerevisiae with well positioned poris, we have checked the robustness of our results with respect to a stochastic number of firing factors ${N}_{D}^{T}$ from cell to cell (Poisson distribution, IyerBiswas et al. (2009)). We confirmed the $I\left(t\right)$ bell shape with a robust duration of the Sphase of $58.6\pm 4.3$ min as compared to $58.5\pm 3.3$ min obtained previously with a constant number of firing factors. Interestingly, in an experiment where hydroxyurea (HU) was added to the yeast growth media, the sequence of activation of replication origins was shown to be conserved even though ${T}_{phase}$ was lengthened from 1 hr to 16 hr (Alvino et al., 2007). HU slows down the DNA synthesis to a rate of $\sim 50\text{\hspace{0.17em}bp\hspace{0.17em}}{\text{min}}^{1}$ corresponding to a 30fold decrease of the fork speed (Sogo et al., 2002). Up to a rescaling of time, the replication kinetics of our model is governed by the ratio between replication fork speed and the productiveinteraction rate ${k}_{on}$ (neglecting here the possible contribution of the activation dynamics of firing factors). Hence, our model can capture the observation of Alvino et al. (2007) when considering a concomitant fork slowing down and ${k}_{on}$ reduction in response to HU, which is consistent with the molecular action of the replication checkpoint induced by HU (Zegerman and Diffley, 2010). In a model where the increase of $I\left(t\right)$ results from the import of replication factors, the import rate would need to be reduced by the presence of HU in proportion with the lengthening of Sphase in order to maintain the pattern of origin activations. Extracting $I\left(t\right)$ from experimental replication data for cells grown in absence (HU${}^{}$) or presence (HU${}^{+}$) (Alvino et al., 2007), we estimated ${I}_{max}^{\text{HU}}\sim 6.0\text{\hspace{0.17em}}{\text{Mb}}^{1}{\text{min}}^{1}$ and ${I}_{max}^{\text{HU}+}\sim 0.24\text{\hspace{0.17em}}{\text{Mb}}^{1}{\text{min}}^{1}$ for HU${}^{}$ and HU${}^{+}$ cells, respectively. The ratio ${I}_{max}^{\text{HU}}/{I}_{max}^{\text{HU}+}\simeq 24.8\sim {v}^{\text{HU}}/{v}^{\text{HU}+}$ is quite consistent with the prediction of the scaling law (Equation (5)) for a constant density of poris.
D. melanogaster and human
We gathered from the literature experimental estimates of ${I}_{max}$, ${\rho}_{0}$ and $v$ for different eukaryotic organisms (Table 1). As shown in Figure 2 (c), when plotting ${I}_{max}$ vs $v{\rho}_{0}^{2}$, all the experimental data points remarkably follow the diagonal trend indicating the validity of the scaling law (Eq. (5)) for all considered eukaryotes. We performed two series of simulations for fixed values of parameters ${k}_{o}$, ${N}_{D}^{T}$ and $v$ and decreasing values of ${\rho}_{0}$ with both periodic distribution (blue) and uniform (green) distributions of poris (Figure 2 (c)). The first set of parameters was chosen to cover high ${I}_{max}$ values similar the one observed for Xenopus embryo (bullets, solid lines). When decreasing ${\rho}_{0}$, the number of firing factors becomes too large and $I\left(t\right)$ does no longer present a maximum. We thus decreased the value of ${N}_{D}^{T}$ keeping all other parameters constant (boxes, dashed line) to explore smaller values of ${I}_{max}$ in the range of those observed for human and D. melanogaster. We can observe that experimental data points’ deviation from Equation (5) is smaller than the deviation due to specific poris distributions.
Note that in human it was suggested that early and late replicating domains could be modeled by spatial inhomogeneity of the pori distribution along chromosomes, with a high density in early replicating domains (${\rho}_{0,early}=2.6$ ORC/100 kb) and a low density in late replicating domains (${\rho}_{0,late}=0.2$ ORC/100 kb) (Miotto et al., 2016). If low and highdensity regions each cover one half of the genome and ${\rho}_{0,early}\gg {\rho}_{0,late}$, most poris are located in the highdensity regions and the origin firing kinetics (${N}_{fired}(t,t+dt)$) will mainly come from initiation in these regions. However, the length of unreplicated DNA also encompasses the late replicating domains resulting in a lowering of the global $I\left(t\right)$ by at least a factor of 2 (Equation (1)). Hence, in the context of our model ${I}_{max}\lesssim 0.5v{\rho}_{early}^{2}$. Interestingly, considering the experimental values for the human genome ($I}_{max}=0.3{\text{Mb}}^{1}{\text{min}}^{1$ and $v=1.46{\text{kb min}}^{1}$, Table 1), this leads to ${\rho}_{0,early}\gtrsim 2.3$ Ori/100 kb, in good agreement with the estimated density of 2.6 ORC/100 kb (Miotto et al., 2016). Inhomogeneities in origin density could create inhomogeneities in firing factor concentration that would further enhance the replication kinetics in high density regions, possibly corresponding to early replication foci.
Discussion
To summarize, we have shown that within the framework of 1D nucleation and growth models of DNA replication kinetics (Herrick et al., 2002; Jun and Bechhoefer, 2005), the sufficient conditions to obtain a universal bell shaped $I\left(t\right)$ as observed in eukaryotes are a strong bimolecular reaction between localized poris and limiting origin firing factors that travel with replication forks and are released at termination. Under these conditions, the density of poris naturally decreases by the end of the Sphase and so does ${I}_{S}\left(t\right)$. Previous models in Xenopus embryo (Goldar et al., 2008; Gauthier and Bechhoefer, 2009) assumed that all sites contained a pori implying that the time ${t}_{r}$ to replicate DNA between two neighboring poris was close to zero. This clarifies why they needed some additional mechanisms to explain the final decrease of the firing rate. Moreover, our model predicts that the maximum value for $I\left(t\right)$ is intimately related to the density of poris and the fork speed (Equation (5)), and we have shown that without free parameter, this relationship holds for five species with up to a 300fold difference of ${I}_{max}$ and $v{\rho}_{0}^{2}$ (Table 1, Figure 2 (c)).
Our model assumes that all poris are governed by the same rule of initiation resulting from physicochemically realistic particulars of their interaction with limiting replication firing factors. Any spatial inhomogeneity in the firing rate $I(x,t)$ along the genomic coordinate in our simulations thus reflects the inhomogeneity in the distribution of the potential origins in the genome. In yeast, replication kinetics along chromosomes were robustly reproduced in simulations where each replication origin fires following a stochastic law with parameters that change from origin to origin (Yang et al., 2010). Interestingly, this heterogeneity between origins is captured by the MultipleInitiator Model where origin firing time distribution is modeled by the number of MCM27 complexes bound at the origin (Yang et al., 2010; Das et al., 2015). In human, early and late replicating domains could be modeled by the spatial heterogeneity of the origin recognition complex (ORC) distribution (Miotto et al., 2016). In these models, MCM27 and ORC have the same status as our poris, they are potential origins with identical firing properties. Our results show that the universal bellshaped temporal rate of replication origin firing can be explained irrespective of speciesspecific spatial heterogeneity in origin strength. Note, however, that current successful modeling of the chromosome organization of DNA replication timing relies on heterogeneities in origins’ strength and spatial distributions (Bechhoefer and Rhind, 2012).
To confirm the simple physical basis of our modeling, we used molecular dynamics rules as previously developed for S. cerevisiae (Arbona et al., 2017) to simulate Sphase dynamics of chromosomes confined in a spherical nucleus. We added firing factors that are free to diffuse in the covolume left by the chain and that can bind to proximal poris to initiate replication, move along the chromosomes with the replication forks and be released when two fork merges. As shown in Figure 2 (a,b) for Xenopus embryo and S. cerevisiae, results confirmed the physical relevance of our minimal modeling and the validity of its predictions when the 3D diffusion of the firing factors is explicitly taken into account. Modeling of replication timing profiles in human was recently successfully achieved in a model with both inhibition of origin firing 55 kb around active forks, and an enhanced firing rate further away up to a few 100 kb (Löb et al., 2016) as well as in models that do not consider any inhibition nor enhanced firing rate due to fork progression (Gindin et al., 2014; Miotto et al., 2016). These works illustrate that untangling spatiotemporal correlations in replication kinetics is challenging. 3D modeling opens new perspectives for understanding the contribution of firing factor transport to the correlations between firing events along chromosomes. For example in S. cerevisiae (Knott et al., 2012) and in S. pombe (Kaykov and Nurse, 2015), a higher firing rate has been reported near origins that have just fired (but see Yang et al. (2010)). In mammals, megabase chromosomal regions of synchronous firing were first observed a long time ago (Huberman and Riggs, 1968; Hyrien, 2016b) and the projection of the replication program on 3D models of chromosome architecture was shown to reproduce the observed Sphase dynamics of replication foci (Löb et al., 2016). Recently, profiling of replication fork directionality obtained by Okazaki fragment sequencing have suggested that early firing origins located at the border of Topologically Associating Domains (TADs) trigger a cascade of secondary initiation events propagating through the TAD (Petryk et al., 2016). Early and late replicating domains were associated with nuclear compartments of open and closed chromatin (Ryba et al., 2010; Boulos et al., 2015; Goldar et al., 2016; Hyrien, 2016b). In human, replication timing Udomains (0.1–3 Mb) were shown to correlate with chromosome structural domains (Baker et al., 2012; Moindrot et al., 2012; Pope et al., 2014) and chromatin loops (Boulos et al., 2013, Boulos et al., 2014).
Understanding to which extent spatiotemporal correlations of the replication program can be explained by the diffusion of firing factors in the tertiary chromatin structure specific to each eukaryotic organism is a challenging issue for future work.
Materials and methods
Wellmixed model simulations
Request a detailed protocolEach model simulation allows the reconstruction of the full replication kinetics during one Sphase. Chromosome initial replication state is described by the distribution of poris along each chromosomes. For Xenopus embryo, pori positions are randomly determined at the beginning of each simulation following two possible scenarios:
For the uniform distribution scenario, $L{\rho}_{0}$ origins are randomly positions in the segment $[0,L]$, where ${\rho}_{0}$ is the average density of potential origins and $L$ the total length of DNA.
For the periodic distribution scenario, exactly one origin is positioned in every nonoverlapping $1/{\rho}_{0}$ long segment. Within each segment, the position of the origin is chosen randomly in order to avoid spurious synchronization effects.
For yeast, the pori positions are identical in each Sphase simulations and correspond to experimentally determined positions reported in OriDB (Siow et al., 2012). The simulation starts with a fixed number ${N}_{D}^{T}$ of firing factors that are progressively made available as described in Results. At every time step $t=ndt$, each free firing factor (available factors not bound to an active replication fork) has a probability to fire one of the ${N}_{pori}\left(t\right)$ poris at unreplicated loci given by:
A random number is generated, and if it is inferior to this probability, an unreplicated pori is chosen at random, two diverging forks are created at this locus and the number of free firing factors decreases by 1. Finally, every fork is propagated by a length $vdt$ resulting in an increase amount of DNA marked as replicated and possibly to the passivation of some poris. If two forks meet they are removed and the number of free firing factors increases by 1. Forks that reach the end of a chromosome are discarded. The numbers of firing events (${N}_{fired}\left(t\right)$), origin passivations, free firing factors (${N}_{FD}\left(t\right)$) and unreplicated poris (${N}_{\text{p}\text{ori}}\left(t\right)$) as well as the length of unreplicated DNA (${L}_{unrepDNA}\left(t\right)$) are recorded allowing the computation of ${I}_{S}\left(t\right)$ (Eq. (1)), the normalized density of poris (${\rho}_{\text{p}\text{ori}}\left(t\right))/{\rho}_{0}$), the normalized number of free firing factors (${N}_{FD}\left(t\right)/{N}_{FD}^{\ast}\left(t\right)$) and the ratio between the number of origin passivations and activations. Simulation ends when all DNA has been replicated, which define the replication time.
3D model simulations
Request a detailed protocolReplication kinetics simulation for the 3D model follows the same steps as in the wellmixed model except that the probability that a free firing factor activates an unreplicated pori depends on their distance $d$ obtained from a molecular dynamic simulation performed in parallel to the replication kinetics simulation. We used HOOMDblue (Anderson et al., 2008) with parameters similar to the ones previously considered in Arbona et al. (2017) to simulate chromosome conformation dynamics and free firing factor diffusion within a spherical nucleus of volume ${V}_{N}$. The details of the interaction between the diffusing firing factors and the poris is illustrated in Figure 2—figure supplement 1. Given a capture radius ${r}_{c}$ set to two coarse grained chromatin monomer radiuses, when a free firing factor is within the capture volume ${V}_{c}={\scriptscriptstyle \frac{4}{3}}\pi {r}_{c}^{3}$ around an unreplicated pori ($d<{r}_{c}$), it can activate the origin with a probability $p$. In order to have a similar firing activity as in the wellmixed model, ${r}_{c}$ and $p$ were chosen so that $p{V}_{c}/{V}_{N}$ takes a value comparable to the ${k}_{on}$ values used in the wellmixed simulations.
For each set of parameters of the wellmixed and 3D models, we reported the mean curves obtained over a number of independent simulations large enough so that the noisy fluctuations of the mean ${I}_{S}\left(t\right)$ are small compared to the average bellshaped curve. The complete set of parameters for each simulation series is provided in Supplementary file 1. The scripts used to extract yeast $I\left(t\right)$ from the experimental data of Alvino et al. (2007) can be found here https://github.com/ jeammimi/ifromprof/blob/master/notebooks/exploratory/Alvino_WT.ipynb (yeast in normal growth conditions) and here https://github.com/jeammimi/ifromprof/blob/master/notebooks/exploratory/Alvino_H.ipynb (yeast grown grown in Hydroxyurea) (Arbona and Goldar, 2018). A copy is archived at https://github.com/elifesciencespublications/ifromprof.
Data availability
All experimental data analyzed in this study are included in the manuscript. Source data files have been provided for Figure 2.
References

Replication in hydroxyurea: it's a matter of timeMolecular and Cellular Biology 27:6396–6406.https://doi.org/10.1128/MCB.0071907

General purpose molecular dynamics simulations fully implemented on graphics processing unitsJournal of Computational Physics 227:5342–5359.https://doi.org/10.1016/j.jcp.2008.01.047

Elucidating the DDKdependent step in replication initiationThe EMBO Journal 35:907–908.https://doi.org/10.15252/embj.201694227

Replication timing and its emergence from stochastic processesTrends in Genetics 28:374–381.https://doi.org/10.1016/j.tig.2012.03.011

Structural organization of human replication timing domainsFEBS Letters 589:2944–2957.https://doi.org/10.1016/j.febslet.2015.04.015

Replication fork velocities at adjacent replication origins are coordinately modified during DNA replication in human cellsMolecular Biology of the Cell 18:3059–3067.https://doi.org/10.1091/mbc.e06080689

DNA combing reveals intrinsic temporal disorder in the replication of yeast chromosome VIJournal of Molecular Biology 375:12–19.https://doi.org/10.1016/j.jmb.2007.10.046

Replication timing is regulated by the number of MCMs loaded at originsGenome Research 25:1886–1892.https://doi.org/10.1101/gr.195305.115

Mathematical modelling of whole chromosome replicationNucleic Acids Research 38:5623–5633.https://doi.org/10.1093/nar/gkq343

Control of DNA replication by anomalous reactiondiffusion kineticsPhysical Review Letters 102:158104.https://doi.org/10.1103/PhysRevLett.102.158104

A chromatin structurebased model accurately predicts DNA replication timing in human cellsMolecular Systems Biology 10:722.https://doi.org/10.1002/msb.134859

Kinetic model of DNA replication in eukaryotic organismsJournal of Molecular Biology 320:741–750.https://doi.org/10.1016/S00222836(02)005223

On the mechanism of DNA replication in mammalian chromosomesJournal of Molecular Biology 32:327–341.https://doi.org/10.1016/00222836(68)900132

Mathematical modelling of eukaryotic DNA replicationChromosome Research 18:147–161.https://doi.org/10.1007/s1057700990924

From simple bacterial and archaeal replicons to replication N/UdomainsJournal of Molecular Biology 425:4673–4689.https://doi.org/10.1016/j.jmb.2013.09.021

BookUp and down the slope: Replication timing and fork directionality gradients in eukaryotic genomesIn: Kaplan D. L, editors. The Initiation of DNA Replication in Eukaryotes. Switzerland: Springer International Publishing. pp. 65–85.

Stochasticity of gene products from transcriptional pulsingPhysical Review E 79:031911.https://doi.org/10.1103/PhysRevE.79.031911

Cell cycle regulation of the replication licensing system: involvement of a Cdkdependent inhibitorThe Journal of Cell Biology 136:125–135.https://doi.org/10.1083/jcb.136.1.125

3D chromatin conformation correlates with replication timing and is conserved in resting cellsNucleic Acids Research 40:9470–9481.https://doi.org/10.1093/nar/gks736

Replication landscape of the human genomeNature Communications 7:10208.https://doi.org/10.1038/ncomms10208

OriDB, the DNA replication origin database updated and extendedNucleic Acids Research 40:D682–D686.https://doi.org/10.1093/nar/gkr1091
Article and author information
Author details
Funding
Institut National Du Cancer (PLBIO16302)
 Olivier Hyrien
 Benjamin Audit
Fondation pour la Recherche Médicale (DEI20151234404)
 Arach Goldar
 Olivier Hyrien
 Benjamin Audit
Agence Nationale de la Recherche (ANR15CE12001101)
 Olivier Hyrien
 Alain Arneodo
 Benjamin Audit
Joint Research Institute for Science and Society (JoRISS 20172018)
 Benjamin Audit
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank F. Argoul for helpful discussions. This work was supported by Institut National du Cancer (PLBIO16302), Fondation pour la Recherche Médicale (DEI20151234404) and Agence Nationale de la Recherche (ANR15CE12001101). BA acknowledges support from Science and Technology Commission of Shanghai Municipality (15520711500) and Joint Research Institute for Science and Society (JoRISS). We gratefully acknowledge support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. We thank BioSyL Federation and Ecofect LabEx (ANR11LABX0048) for inspiring scientific events.
Version history
 Received: January 18, 2018
 Accepted: May 31, 2018
 Accepted Manuscript published: June 1, 2018 (version 1)
 Version of Record published: July 5, 2018 (version 2)
Copyright
© 2018, Arbona et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,827
 Page views

 291
 Downloads

 8
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Chromosomes and Gene Expression
 Genetics and Genomics
Heterogeneity in endothelial cell (EC) subphenotypes is becoming increasingly appreciated in atherosclerosis progression. Still, studies quantifying EC heterogeneity across whole transcriptomes and epigenomes in both in vitro and in vivo models are lacking. Multiomic profiling concurrently measuring transcriptomes and accessible chromatin in the same single cells was performed on six distinct primary cultures of human aortic ECs (HAECs) exposed to activating environments characteristic of the atherosclerotic microenvironment in vitro. Metaanalysis of singlecell transcriptomes across 17 human ex vivo arterial specimens was performed and two computational approaches quantitatively evaluated the similarity in molecular profiles between heterogeneous in vitro and ex vivo cell profiles. HAEC cultures were reproducibly populated by four major clusters with distinct pathway enrichment profiles and modest heterogeneous responses: EC1angiogenic, EC2proliferative, EC3activated/mesenchymallike, and EC4mesenchymal. Quantitative comparisons between in vitro and ex vivo transcriptomes confirmed EC1 and EC2 as most canonically EClike, and EC4 as most mesenchymal with minimal effects elicited by siERG and IL1B. Lastly, accessible chromatin regions unique to EC2 and EC4 were most enriched for coronary artery disease (CAD)associated singlenucleotide polymorphisms from Genome Wide Association Studies (GWAS), suggesting that these cell phenotypes harbor CADmodulating mechanisms. Primary EC cultures contain markedly heterogeneous cell subtypes defined by their molecular profiles. Surprisingly, the perturbations used here only modestly shifted cells between subpopulations, suggesting relatively stable molecular phenotypes in culture. Identifying consistently heterogeneous EC subpopulations between in vitro and ex vivo models should pave the way for improving in vitro systems while enabling the mechanisms governing heterogeneous cell state decisions.

 Chromosomes and Gene Expression
We have developed a deep sequencingbased approach, RecSeq, that allows simultaneous monitoring of ribosomal 48S preinitiation complex (PIC) formation on every mRNA in the translatome in an in vitro reconstituted system. RecSeq isolates key early steps in translation initiation in the absence of all other cellular components and processes. Using this approach, we show that the DEADbox ATPase Ded1 promotes 48S PIC formation on the start codons of >1000 native mRNAs, most of which have long, structured 5′untranslated regions (5′UTRs). Remarkably, initiation measured in RecSeq was enhanced by Ded1 for most mRNAs previously shown to be highly Ded1dependent by ribosome profiling of ded1 mutants in vivo, demonstrating that the core translation functions of the factor are recapitulated in the purified system. Our data do not support a model in which Ded1acts by reducing initiation at alternative start codons in 5′UTRs and instead indicate it functions by directly promoting mRNA recruitment to the 43S PIC and scanning to locate the main start codon. We also provide evidence that eIF4A, another essential DEADbox initiation factor, is required for efficient PIC assembly on almost all mRNAs, regardless of their structural complexity, in contrast to the preferential stimulation by Ded1 of initiation on mRNAs with long, structured 5′UTRs.