Cell-to-cell variability and robustness in S-phase duration from genome replication kinetics
Qing Zhang Federico Bassetti Marco Gherardi Marco Cosentino Lagomarsino


APS/123-QED





 
  qzhang519@gmail.com

 federico.bassetti@unipv.it

 gocram@gmail.com

 marco.cosentino-lagomarsino@upmc.fr

 ^1 Sorbonne
  Universités, UPMC Univ Paris 06, UMR 7238, Computational and
  Quantitative Biology, 15 rue de l'École de Médecine Paris,
  France
and
^2Dipartimento di Matematica, Università di Pavia, Pavia, Italy
^3IFOM, FIRC Institute of Molecular Oncology, Milan, Italy
^4CNRS, UMR 7238, Paris, France



             


  Genome replication, a key process for a cell, relies on stochastic
  initiation by replication origins, causing a variability of
  replication timing from cell to cell.  While stochastic models of
  eukaryotic replication are widely available, the link between the
  key parameters and overall replication timing has not been addressed
  systematically.

  We use a combined analytical and computational approach to calculate how
  positions and strength of many origins lead to a given cell-to-cell
  variability of total duration of the replication of a large region,
  a chromosome or the entire genome.

  Specifically, the total replication timing can be framed as an
  extreme-value problem, since it is due to the last region that
  replicates in each cell.

  Our calculations identify two regimes based on the spread between
  characteristic completion times of all inter-origin regions of a
  genome. For widely different completion times, timing is set by the
  single specific region that is typically the last to replicate in
  all cells.

  Conversely, when the completion time of all regions are comparable,
  an extreme-value estimate shows that the cell-to-cell variability of
  genome replication timing has universal properties. Comparison with
  available data shows that the replication program of three yeast
  species falls in this extreme-value regime.



                             

                              
Cell-to-cell variability and robustness in S-phase duration from
  genome replication kinetics 
    Marco Cosentino Lagomarsino ^1,3,4
    December 30, 2023
===============================================================================================







§ INTRODUCTION



In all living systems, the duration of DNA replication correlates with
key cell-cycle features, and is intimately linked with transcription,
chromatin structure and genome evolution. Dysfunctional replication
kinetics is associated to cancer and found in aging cells.

Eukaryotic organisms rely on multiple discrete origins of replication
along the DNA <cit.>. These origins
are “licensed” during the G1 phase by origin recognition complexes
and MCM helicases, and can initiate replication during S
phase <cit.>. Once one origin is activated
(“fires”), a pair of replication forks are assembled and move
bidirectionally. In one cell cycle, one origin already activated or
passively replicated cannot be activated
again <cit.>.

Origins have specific firing rates, possibly connected to the
number of bound MCM helicase complexes <cit.>, and their
specificity determines the kinetics of replication during S phase, or
“replication program”.


To investigate genomic replication kinetics, DNA copy number can be
measured with microarray or sequencing, as a function of genome
position and time (see, e.g.,
<cit.>). Based
on such high-throughput replication timing data, it is possible to
infer origin positions and the key parameters for a mathematical
description of the replication process (see,
e.g., <cit.>). Recent
methods also allow to extract the same information from free-cycling
cells <cit.>.  The mathematical modeling of genome-wide
replication timing data shows that replication kinetics results from
the stochastic mechanism of origin
firing <cit.>.

In other words, replication timing originates from individual
probabilities of origin firing (and their correlations with genome
state <cit.>). In such
models, firing rate of individual origins determine the kinetic
pattern of replication along the chromosomal coordinate, and fork
velocity is typically assumed to be nearly constant along the genome
(in absence of blockage).


Evidence of this stochasticity directly from single cells (which
should give access to relevant correlation patterns) is less abundant.
Importantly, replication timing patterns observed in population
studies can be explained by stochastic origin firing at the
single-cell level <cit.>.

Stochastic activation of origins leads to stochasticity of termination
and cell-to-cell variability of the total duration of replication of a
chromosome, a genomic region, or the whole
S-phase <cit.>, with possible repercussions on the
cell cycle. This raises several questions, including how the
individual rates and spatial distribution of origins cooperate to generate
variability in replication timing, the extent of such variability, and
whether it is possible to identify specific regimes or optimization
principles in terms of cell-to-cell variability.

However, such questions have not been systematically addressed in the
available models.




-65.1pt



A series of pioneering studies <cit.> has used
techniques of extreme-value theory to derive the distribution of
replication times in the particular case where each locus of the

genome is a potential origin of replication, as in the embryonic cells
of X. laevis.  These efforts allowed to clarify the possible
optimization principles underlying the replication kinetics in such
organisms.


Here, we extend this approach to the widely relevant case of discrete
origins with fixed
positions <cit.> using
a modeling framework for stochastic replication to investigate the
cell-to-cell variability of the duration of S-phase (or of the
replication of any genomic region such as one chromosome).  We

use analytical calculations based on extreme-value theory and
simulations,

employ experimental data to infer replication parameters and identify
the main features of empirical origin strengths and positions, and
their response to specific changes.




§ MATERIALS AND METHODS




 §.§ Model


 We make use of a one-dimensional nucleation-growth
model <cit.> of stochastic replication kinetics
with discrete origin locations x_i, similar to models available in
the literature <cit.>.

Activation of origins (firing) is stochastic, and is described as a
non-stationary Poisson process.  The firing rate A_i(t) of the
origin located at x_i is a function of time, A_i(t)=λ_i
t^γθ(t), where θ(t) is the step function, and
λ_i and γ are
constants <cit.>.
We assume that the parameter γ and the fork velocity v are
common to all origins, whereas λ_i, which reflects the
specific strength of each origin, is origin dependent.  The
probability density function (PDF) f_i(t) of the firing time t for
the i-th origin, given that the origin fires during that
  replication round, can be obtained as f_i(t)=A_i(t)
exp(-∫_0^t A_i(τ) dτ), which gives

    f_i(t)=λ_i   t^γ θ(t)
           exp(-λ_i t^γ+1/γ+1).

When γ>0, i.e., when the firing rate increases with time,
f_i(t) is a stretched exponential distribution.

When γ=0, the firing rates are constant and the process is
stationary, so A_i(t)=λ_i and f_i(t)=λ_i
θ(t)e^-λ_i t.

Once an origin has fired, replication forks proceed bidirectionally at
constant speed, possibly overriding other origins by passive
replication.  When two forks meet in an inter-origin region,
replication of that region is terminated.

The length of the i-th region is defined as d_i=x_i+1-x_i;
the time when its replication is completed is T_i.

The duration of the S phase T_S
is the time needed for all inter-origin regions to be replicated.






 §.§ Fits


Empirical parameters were inferred through fitting experimental data
from refs. <cit.> on DNA copy
number as a function of position and time with the model.  The
positions of replication origins were obtained directly from the
literature and considered
fixed <cit.>. The fits are
performed by minimizing the distance between the replication timing
profiles in the model and in the experimental data. This is carried
out by updating the global parameters (γ and v) and the local
parameters (λ_i, i∈{1,2,...,n}) iteratively (Appendix A). The parameters from these fits are presented
in Supplementary Table S1.




 §.§ Simulations


Our theoretical calculations (described below)
allow to obtain the cell-to-cell variability of T_S in
special regimes.


We compare simulations using the complete information on the locations
and strengths of all origins fitted from the data, with randomized
chromosomes having similar properties. In these randomized chromosomes
we consider the inter-origin distances d_i and the strengths
λ_i as independent random variables.  They are drawn from
probability distributions recapitulating their empirical mean and
variability.

More precisely, from the fitted parameters we fix the mean <
  d> and the standard deviation σ_d of the distance, and
the mean <λ> and the standard deviation
σ_λ of the strength.

The actual distances d_i and strengths λ_i are then drawn by
sampling from two gamma distributions

    d_i∼Γ(<d>^2/σ_d^2,<d>/σ_d^2),
        λ_i∼Γ(<λ>^2/σ_λ^2,<λ>/σ_λ^2).



The gamma distribution Γ(a,b)
(parametrized in terms of a shape parameter a and a rate parameter b)
has PDF p(x) ∝ x^a-1exp(-b x).  It yields positive values, with mean a/b and variance
a/b^2, and it is the maximum-entropy distribution with fixed mean
and fixed mean of the logarithm. We verified that the assumption of a
gamma distribution was in line with empirical data (Fig. S1).



To explore the full range of parameters, we also used stochastic
simulations, which were performed both (i) with the precise origin
locations and strengths fitted from the data, and (ii) with d_i and
λ_i drawn randomly as described above.


To avoid the boundary effects of linear chromosomes, we consider
circular chromosomes with n origins, unless specified otherwise
(boundary effects are discussed in the Appendix B and Fig. S2,
and do not affect our main conclusions.)



To analyze the biologically relevant regimes, we considered
replication kinetics data on different yeast species, from
refs. <cit.> and
    <cit.>, ran simulations with such
parameters, and compared with the theoretical predictions using the
empirical values for σ_d, σ_λ and mean origin
positions and strengths.




§ BACKGROUND







 §.§ The S-phase duration is the result of a maximum
operation on the stochastic replication times of inter-origin regions



We start by discussing how the stochastic nature of single-origin firing
affects the total replication timing of a chromosome.
Fig. <ref>ab illustrates this process.  In each cell, a
chromosome is fully replicated when the last inter-origin region is
complete.  In other words, the last-replicated region sets the
completion time for the whole chromosome.  Consequently, the total
duration is the maximum among the replication times of all
inter-origin regions <cit.>.

For simplicity, we first consider the case of a genome with only one
chromosome.

The duration of the S phase is therefore
T_S=max(T_1,T_2,...,T_n) where n is the number of inter-origin
regions.

The stochasticity of the replication time T_i of each inter-origin
region makes the S-phase duration T_S itself stochastic,
thus giving rise to cell-to-cell variability, which can be estimated
by the model (Fig. <ref>c). In the case of multiple
chromosomes, the same reasoning applies to the last-replicated
inter-origin region over all chromosomes.




 






§ RESULTS





 §.§ A theoretical calculation reveals the existence of two
  distinct regimes for the replication program


  

It is possible to estimate the distribution of T_S
analytically, starting from the distribution of T_i.

Two distinct limit-case scenarios can be distinguished.  In the first
scenario, a specific inter-origin region r is typically the slowest
to complete replication and thus represents a “replication
bottleneck”.  In this case, T_S is dominated by T_r,
meaning that T_S≈ T_r.

T_r is identified as the one which is largest on average.

Fig. <ref>a shows an example chromosome with 10 origins
with the same strength, where one inter-origin distance (d_1) is
much larger than the others.

Owing to this disparity, T_1 is very likely  the maximum among
all T_i, and is therefore the region determining T_S.

In this scenario, which we term “bottleneck estimate”, the
distribution of T_S will be approximately the same as that
of the bottleneck T_r (Fig. <ref>c).



In the second scenario, each inter-origin region has a similar
probability to be the latest to complete replication.

In this case, every inter-origin region contributes to the distribution
of T_S.  Since T_S=max(T_1,T_2,…,T_n), we apply the
well-known Fisher-Tippett-Gnedenko theorem
<cit.>, which is a general
result on extreme-value distributions (EVD).

In order to use this theorem, we make the following two assumptions:
(i) T_1, T_2, …, T_n are statistically independent, i.e.,
each inter-origin replication time is an independent random variable,
incorporating the essential information about origin variability and rates;
(ii) T_i follows a stretched-exponential distribution,
independent of i, i.e.

    p(T_i<t)=1-e^-α(t-t_0)^β,

when t>t_0, while p(T_i<t)=0 when t⩽ t_0. The (positive)
parameters α, β and t_0, effectively describe the
consequences of the model parameters v, γ, inter-origin
distances (d_1, d_2, ..., d_n) and origin strengths
(λ_1,λ_2,...,λ_n) on completion timing of
inter-origin regions (see below and Appendix D), and can be
obtained by fitting the distribution of replication time for a typical
inter-origin region (obtained from simulations) with Eq. <ref>.


Our fits show that Eq. <ref> is a remarkably good phenomenological
approximation of the distribution of T_i (see Appendix C and
Fig. S3), thus justifying assumption (ii) above.

Note that the fitted stretched exponential form also incorporates
effectively the coupling existing between different inter-origin
regions.

Indeed, neighboring regions are correlated since they use a pair of
replication forks stemming from their common origin.  Moreover, even
distant inter-origin regions can share the same fork if they are
passively replicated.

In order to justify the assumption (i), we tested the effect of the
correlation between different regions, by sampling T_1, T_2, …,
T_n from the distribution in Eq. <ref> independently and then
taking their maximum T_S^*. We verified that the difference
between the distribution of T_S^* and that of
T_S obtained from simulation (where the correlations are
present) is small.  Therefore, the effect of these relatively
  short-ranged correlations can be, to a first approximation,
  neglected at the scale of the chromosomes and of the genome, and
  described by the effective stretched-exponential form (see
Fig. S4).


Based on these assumptions, we can use the Fisher-Tippett-Gnedenko
theorem and derive the following cumulative distribution function for
T_S as a function of the number of origins n
and the parameters α, β and t_0 (the
calculation is detailed in the Appendix D):

    P(T_S≤ t)≈exp{-exp[βlog
          n(1-(α/log n)^1/β(t-t_0))]}.



 Eq. <ref> gives a direct estimate of
the distribution of the S-phase duration in this second scenario,
which we term “extreme-value” or “EVD” regime.
The resulting
distribution is universal, since it does not depend on the detailed
positions and rates of the origins, and depends in a simple way on the
parameters α, β, t_0 and n.

Although the extreme-value estimate should apply to the case of
  large n, the approximation Eq. <ref> holds to a satisfactory
  extent also for realistic values, when n is order 10 (see
  Supplementary Fig. S12).

We also derived approximate analytical expressions for α,
β and t_0 as functions of the parameters v, γ, for a
“typical” region characterized by <λ> and
<d> under the assumption of negligible interference from
non-neighbour origins (see Appendix D).



The procedure by which we apply Eqs. <ref> and <ref> is the
following. Given inter-origin distances and origins strengths assigned
arbitrarily or inferred from empirical data, the simulation of the
replication of a chromosome gives the distribution of T_i and
T_S. A fit of the distribution of T_i from simulation
using Eq. <ref> gives the parameters α, β and t_0. Finally,
the EVD estimate for the distribution of T_S, can be
obtained from Eq. <ref>, and compared with the distribution of
T_S form simulations.

This procedure can be seen as a variant of the method introduced
in refs. <cit.> applicable to the case of discrete
origins (see Discussion).




Fig. <ref>b shows one example where one circular chromosome
has 10 origins with identical strengths and identical inter-origin
distances. The estimated distribution of S-phase duration from
Eq. <ref> is well-matched with the simulated one (Fig. <ref>d).

Fig. <ref> also shows how the bottleneck estimate works for
the opposite scenario, and compares simulations with both estimates in
the two different regimes.

Similar to Fig. <ref>, Supplementary Fig. S5 shows the
existence of the two regimes in presence of a single origin affecting
the two neighboring inter-origin regions.  In the bottleneck regime,
these two regions replicate much later than the others, because their
common origin is much weaker than the other origins; the S-phase
duration is then dominated by their replication time. This case also
illustrates how the bottleneck regime may not be limited to a single
inter-origin region.

Finally, Supplementary Fig. S6 shows the distribution of the
inter-origin completion times T_i in the cases presented in
Fig. <ref> and Supplementary Fig. S5. This analysis
illustrates how extra peaks in the right tail of T_i distribution
relate to the failure of the extreme-value estimate for the
distribution of S-phase duration. These examples indicate that, as
expected, the presence of outliers in the values of T_i
(exceedingly slowly-replicating regions) is responsible for the
onset of the bottleneck behavior.





 §.§ The extreme-value regime is robust to perturbations
  increasing the replication timing of a local region


Origin number, origin strengths and inter-origin distances can be
perturbed due to genetic change (DNA mutation or recombination), over
evolution, and due to epigenetic effects such as binding of specific
agents.  We can compare the robustness of the two regimes identified
above to perturbations of these parameters.  We consider in particular
the elongation of a single inter-origin distance d_i ↦ d_i +
δ_d (similar results to those reported below are obtained
  for a perturbation affecting the strength of a single origin, see
  Supplementary Fig. S7).


In such case, the change of T_i is approximately equal to δ_d/2v.

In the bottleneck regime, if the perturbed inter-origin region is
the slowest-replicating one, <T_S>
increases linearly with δ_d with slope 1/2v, and the
distribution of T_S shifts by a delay δ_d/2v
(Fig. <ref>a).

In the extreme-value regime, instead, there is no single bottleneck
inter-origin region, and the change of
T_S with the perturbation turns out to be much smaller than
δ_d/2v (Fig. <ref>b).


Notice that in both regimes the variability of the S-phase duration
around its average is not affected sensibly (insets of
Fig. <ref>).

In summary, the bottleneck

region is “sensitive”

to the specific perturbations considered, since termination of
replication is highly dependent on a single inter-origin region, while
the EVD regime is “robust”, as the effect of small local
perturbations can be absorbed by passive replication from nearby
  origins <cit.>.









 §.§ Diversity between completion times of inter-origin
  regions sets the regime of the replication program


The cases discussed above (Fig. <ref>) recapitulate the
expected behavior  in case of high versus small variability of
the typical completion time of different inter-origin regions.

One can expect that if the variability of the inter-origin distances
is large, or origin strengths are heterogenous, it will be more likely
to produce a bottleneck region, which in turn will trivially affect
replication timing.

Conversely, the replication program will be in the extreme-value regime
if the completion times of all regions are comparable.

In order to show this, we tested systematically how average and
variability of T_S change with the variability of
inter-origin distances and origin strengths in randomly generated
genomes.  In this analysis, origin spacings and strengths are
  assigned according to the prescribed probability distributions shown
  in Eq. <ref>, with varying parameters (see the Methods for a
  precise description of how chromosomes are generated).





Fig. <ref> shows the results.  Importantly, we find that
the regimes defined above as extreme cases apply for most
  parameter sets, and there is only a small region of the parameters
  where we find intermediate cases.

Specifically, two parameters, the standard deviations σ_d and
σ_λ, of the inter-origins distances and the origin
strengths respectively, are sufficient to characterize the system.

Fig. <ref>a indicates that as long as σ_d is smaller
than a threshold (around 30 kb), the average
<T_S> and the standard deviation
σ(T_S) of the replication time are approximately
constant.  In this regime, the extreme-value estimate matches well the
simulation results.  When σ_d exceeds the threshold, the
average of T_S increases and its standard deviation
decreases with large fluctuations.  In this other regime, both
<T_S> and σ(T_S) deviate from
the EVD estimate.

Fig. <ref>b shows that varying σ_λ at fixed
origin positions produces a similar behavior (although with smaller
deviations from the EVD estimates).




This analysis shows an emergent dichotomy between these two regimes,
which depends on the distribution of T_i (i.e. both inter-origin
distances and origin firing rates). In principle, more complex
situations where e.g. a subset of many comparably “slow”
inter-origin regions dominates S-phase timing is possible, but this
situation is very rare (and negligible) if origin rates and positions
are generated with the criteria used here (given by
Eq. <ref>).

De facto, under these prescriptions, motivated by empirical
properties of origin positions and strengths, only the two regimes
defined above as extreme cases were observable.  For example, one
  can imagine a situation where each chromosome are, separately, in
  the EVD regime, but the replication of one of the chromosomes takes
  considerably longer than the others on average, which may lead the
  S-phase duration to be in the bottleneck regime. However, we find
  that this situation is essentially never found if origin rates and
  positions have empirically relevant values (i.e. for all
  realizations with empirical means and variances of inter-origin
  distances and origin firing rates).

Qualitatively, this will always be the case if the distribution of
T_i shows a single mode, and there are very few, or just one
  exceptional late-replicating region.



This behavior suggests to define “critical values” of σ_d
and σ_λ, separating the extreme-value regime from the
bottleneck regime, as follows.

We define the σ^c_d, at fixed σ_λ, as the
value of σ_d at which <T_S> (possibly averaged
over many samples of the origin configuration too, denoted
<<T_S>>) is 20% larger than at σ_d=0
and σ_λ=0.  The results presented here do not depend
appreciably on this threshold and do not change much if we define
σ^c_d as the value of σ_d at which
<T_S> is 20% off the prediction of the EVD theory.
The same definition holds for σ^c_λ at fixed
σ_d.

Surprisingly, σ^c_d turns out to be independent of
σ_λ, and σ^c_λ independent of
σ_d.
The resulting “phase diagram”, shown in
Fig. <ref>c, separates the space of parameters into an approximately
rectangular region where the EVD estimate is precise, and an
outer region where heterogeneities dominate, which is identified with
the bottleneck regime.


We can give a simple argument for why this phase diagram is approximately
rectangle-shaped.  Intuitively, a large σ_d increases the
probability of extracting a very large value for d, and a large
σ_λ increases the probability of extracting a very small
λ.  In a realization of a randomized chromosome, such rare
events may generate an extremely slow-replicating region acting as the
bottleneck.

Clearly, drawing an extreme value for only one of the two variables is
sufficient to generate the bottleneck region, giving rise to the two
sides of the rectangle.

For values of the variances of both variables that are below the
individual thresholds, drawing a large d and small λ
jointly makes the upper-right region of the rectangle rounded.
However, such joint extreme draws in the same
inter-origin region are very rare, because the two variables are drawn
independently, so the rounded upper-right corner is very small, as
visible in Fig. <ref>c.

















 §.§ The yeast replication program is just inside the EVD
  regime and likely under selection for short S-phase duration



The results of the previous section indicate that the standard
deviations of the origin distances and of the strengths are the most
relevant parameters determining the regime of the distribution of the
S-phase duration across cells.

We inferred the parameters from replication timing data of the yeasts
S. cerevisiae (ref. <cit.>),
L. kluyveri (ref. <cit.>) and S. pombe (ref. <cit.>). Such fits
fully constrain the model parameters: fork velocity v, γ,
start of the S phase t_0, origin strengths λ_i and
inter-origin distances d_i, from which we calculated
<d>, <λ>, σ_d and
σ_λ, and simulated the duration of S phase and
replication time of each chromosome (see Appendix A and
Fig. S8-10).

In these simulations we consider circular chromosomes with n
  origins, and boundary effects are tested in the Appendix B
  and Fig. S2, and do not affect our main conclusions, indicating
  that, according to the model, the partition of the genome into 16
  unconnected chromosomes has little effect on the statistics of
  S-phase duration.

The values of γ that were obtained as best fits of the
empirical data (Supplementary Fig. S8) were in line with previous
analyses (e.g. <cit.>). In addition,
we found that the standard deviation of the predicted S-phase duration
decreases with the parameter γ (Supplementary Fig. S9), which
agrees with the finding of previous studies focused on
X. laevis <cit.>.



This analysis indicates that the whole-genome values of σ_d and
σ_λ measured for S. cerevisiae, L. kluyveri
and S. pombe place these genomes within the extreme-value
regime.  Rescaling σ_d and σ_λ by the crossover
values σ^c_d and σ^c_λ
respectively makes it possible to compare data with different mean
T_S. This comparison (Fig. <ref>a) shows that not
only the genomic but also most of chromosomal parameters of
L. kluyveri, S. cerevisiae and S. pombe are
located in the extreme-value regime.

With the fitted parameters, most of chromosomes and genomes are
  found in the extreme-value regime (as an example, see Supplementary
  Fig.S10).


Interestingly, all chromosomes (and the full genome) lie close to the
transition line.  This may be a consequence of the presence of
competing optimization goals, such as replication speed (or
reliability) and resource consumption by the replication
machinery <cit.>.

Furthermore, we considered data of two S. cerevisiae
  mutants.  In one mutant, three specific origins in three different
  chromosomes (6, 7, and 10) were
  inactivated <cit.>.

  The inactivation of a specific origin slows down the replication of
  the nearby region, which might cause a bottleneck.

  Our results show that this origin mutant is still in EVD regime
  (Supplementary Fig. S13).

  Importantly, in this case the model should be able to make a precise
  prediction for the replication profile of the chromosomes where one
  origin is inactivated.

  Supplementary Fig. S14 shows the prediction on the replication
  profile of origin mutant strain based on the parameters fitted from
  the data of wild-type strain (except that the three inactivated
  origins are deleted from the origin list). The model prediction is
  in fairly good agreement with data. The mismatch between prediction
  and data in some regions (but not others) is an interesting feature
  revealed by the model, and may result from experimental error or
  gene-expression adaptation of the mutants <cit.>.

  The other mutant strain that we considered is isw2/nhp10,
  from the study of Vincent and coworkers <cit.>, who
  analyzed the functional roles of the Isw2 and Ino80 complexes in DNA
  replication kinetics under stress. This study compares the behavior
  of wild type (wt) strain and a isw2/nhp10 mutant in the
  presence of MMS (DNA alkylating agent methyl methanesulfonate) and
  found that S-phase in isw2/nhp10 is extended compared to the
  wt strain because the Isw2 and Ino80 complexes facilitate
  replication in late-replicating-regions and improve replication fork
  velocity. In agreement with these findings, the model fit of the
  data shows that isw2/nhp10 mutant has more inactive origins
  and smaller fork velocity. Such conditions may facilitate the onset
  of a bottleneck regime in the mutant compared to the wt strain. We
  found that S. cerevisiae wt strain treated with MMS still
  falls in the extreme-value regime. Conversely, some chromosomes (e.g
  13 and 15) of the isw2/nhp10 mutant are in the bottleneck
  regime, and in this case, the whole genome (entire S-phase), is
  driven in the bottleneck regime (see Supplementary Fig. S15).

  Strikingly, the model makes a good prediction on the replication
  profile of the isw2/nhp10 mutant, using origin firing
  strengths and the γ values fitted from the wild-type strain
  experiments, and just adjusting two (global) parameters replication
  speed and an overall factor in all origin firing rates
  (Supplementary Fig. S16). This provides a good cross-validation of
  the applicability of the model in a predictive framework.


A further question is whether we can detect signs of optimization in
the duration of chromosome replication.

Fig. <ref>b compare the S-phase durations obtained from
simulations of the model in two cases: (i) by using the origin
positions and strengths from empirical data (see Supplementary
Fig. S10), and (ii) by using a null model with randomized parameters
(both origin strengths and inter-origin distances) drawn according to
Eq. (<ref>), and preserving the empirical mean and variance.

The results show that for some of the chromosomes the average
replication timing T_S is close to the typical one obtained
from randomized origins (e.g., chromosomes 1,3,5,6,8,11,13 in
S. cerevisiae).  For other chromosomes (e.g., 2,4,7,10,12,15,16 in
S. cerevisiae) the empirical average T_S is instead
very close to the minimum reachable within their ensemble of
randomizations.

Remarkably, chromosomes with higher average replication timing in the
randomized ensemble seem to be more subject to pressure towards
decreasing their average T_S (Supplementary Fig. S11).
This result suggests that the whole replication program may be under
selective pressure for fast replication.










§ DISCUSSION



The core of our results are analytical estimates that capture the
cell-to-cell variability in S-phase duration based on the measurable
parameters of replication kinetics.


Extreme-value statistics has been applied to DNA replication
  before <cit.>, but only to the case of
  organisms like X. laevis, where origin positions are not
  fixed and there is no spatial variability of initiation rates. To
  our knowledge, this method has not been applied systematically to
  fixed-origin organisms such as yeast.





More specifically ref. <cit.> explores the case of a perfect
lattice of equally spaced discrete origins with fixed and equal firing rates,
but does not address the role of the variability of inter-origin
replication times due to randomness in firing rates and inter-origin
distance, which is relevant for fixed-origin organisms.








Another difference is that the authors of
ref. <cit.> derive the coalescence
distribution starting from their model, while here we assume a
stretched-exponential, motivated by data analysis.

Since their distribution is more complex (although the model is
simpler), EVD estimate leads to a formula linking the parameters of
the Gumbel distribution to the initiation parameters in the form of an
implicit equation, that needs to be solved numerically.

Conversely, the assumption that the shape of the distribution of T_i
is given (and estimated from data), gives an explicit relationship
between the parameters describing the T_i distribution and the
Gumbel parameters, leading to simpler formulas and applicability to
the case of discrete origins with different spacings and firing
rates. The parameters of the T_i distribution have then to be
related to the microscopic parameters (See Appendix D).



It is important to note that an approach based on extreme-value
distribution theory is general <cit.>.  Simulations
(including the model used here) are based on specific assumptions that
are often not simple to test and many models on the market use
slightly different assumptions.


Instead, the extreme-value estimates are robust to different shades of
assumptions used in the models available in the literature, and thus
more comprehensive.


Our estimates reveal universal behavior in the distribution of S-phase
duration.  There is a prescribed relation between mean and
variance of S-phase duration, defining a “scaling” behavior for its
distribution. Such universality has been observed in cell-cycle
periods and cell size <cit.>.



Qualitatively, we expect the same universality to hold in a regime
when origins have less than 100% efficiencies, and some may not fire
at all during S-phase.  Origins that fire only in a fraction of the
realizations are accounted for in our simulations, but they entail
second-neighbour effects that are not currently accounted in our
estimates.




There are hundreds of origins in a genome, but our
analysis shows that the relevant parameters to capture the overall
behavior are the means and variances of inter-origin distances and
origin firing rates.

Specifically, we find that two regimes describe most of the
phenomenology, and they depend on the values of these effective
variables.

Importantly, the regimes identified here differ from those
  identified in ref. <cit.>, which just identify a critical
  spacing between discrete (equally spaced) origins, for which
  replication timing starts to be linear with inter-origin distance.


 The notion that the last regions to replicate may tend to be
  different in every cell (our “extreme-value” regime) has been
  proposed already by Hawkins and
  coworkers <cit.>. The opposite regime where some
  specific regions tend to always replicate last ('bottleneck
  region'), has been proposed for mammalian common fragile
  sites <cit.>. Such regions of slow replication,
  pausing and frequent termination have also been described in
  yeast <cit.>.
  These studies make it plausible to think that both extreme-value and
  bottleneck regimes may apply to yeast, despite our analysis based on
  replication kinetics data indicating some pressure towards the
  extreme-value regime.

  Another important case for what concerns replication termination is
  the rDNA locus, which cannot be analyzed in replication kinetics data
  based on microarrays / sequencing data due to its repetitive nature
  (150 identical copies in yeast). However, the large inter-origin
  distances, pseudo-unidirectional replication and epigenetic control
  of origin firing in this locus <cit.> make it a good
  candidate for the last sequence to replicate in yeast. 


Importantly the model used here is similar to a set of previous
  studies, which have tested this approach and validated it with
  experimental data
   <cit.>. Our
  analysis of S-phase duration in single cells is generic, and
  expected to be robust to variations model details.

  The mutant data sets analyzed here also support the predictive power
  of the model in presence of perturbations and parameter changes, and
  hence validate the use of the model in a predictive framework.


Our predictions are compatible with the available values for
  average S-phase duration, which can be roughly estimated through
flow cytometry <cit.>, and
corresponds well to the values obtained by the model (around 60
minutes for S. cerevisiae). Other yeast studies found smaller
values in other conditions <cit.>, which would be
interesting to study with the model.

Additionally, we provide a prediction for the cell-to-cell variability
of S-phase duration, which is an important step of the cell cycle.
Indeed, completion of replication needs to be coordinated with growth
and progression of the cell cycle
stages <cit.>. Cell-to-cell variability in
replication kinetics makes the S phase subject to inherent
stochasticity.

Experimentally, measuring the cell-to-cell variation of the S-phase duration
is a challenge.

While some studies exist using mammalian (cancer) cell
  lines <cit.>, they currently do not have the precision
needed to allow a quantitative match with models.

However, we expect that such measurements will become available in the
near future, thanks to rapidly developing methods of single-cell
biology <cit.>.

Our predictions define some key properties of the replication period
that may be tested with, e.g., single-cell studies in budding yeast,
using the parameters available from replication kinetics studies.  In
this model the S phase is (by itself) a “timer”, so its connection
to cell size homeostasis must be affected by external
mechanisms <cit.>. S-phase duration has been measured
on single E. coli cells, and found to be unlinked to cell
size <cit.>.


Interestingly, our predictions of S-phase duration and
  variability as a function of chromosome copy numbers (Supplementary
  Fig. S12) might apply to cancer cell lines with different levels of
  aneuploidy <cit.>.

Finally, there is the possibility of applying this framework to
describe relevant perturbations <cit.>. This
could also help elucidate how response to DNA damage affects the
replication timing and its variability across cells.


Intriguingly, we also found evidence of bias towards faster
replication in empirical chromosomes compared to randomized ones.

Thus, our overall findings support the hypothesis of a possible
  selective pressure for faster replication, and against bottlenecks.

Other approaches have assumed optimization for faster replication and
looked for optimal origin placement <cit.> or
found other signs of optimality in similar
data <cit.>. Our results are in line with these
findings, and isolate a complementary direction for such optimization.

All these considerations support the biological importance of
replication timing of inter-origin regions and its variability.

However, the sources of the constraints remain an open
question. Clearly, overall replication speed can increase indefinitely
by increasing origin number and initiation rates. However, there are
likely yet-to-be-characterized tradeoffs in these quantities, that
prevent this from happening, and force the system to optimize the
duration of replication in a smaller space of parameters. The
molecular basis for such constraints likely lies at least in part in
the finite resources available for initiation
complexes <cit.>.





We are grateful to Gilles Fischer, Nicolas Agier, Alessandra Carbone
and Renaud Dessalles for useful discussions.

QZ was supperted by the LabEx CALSIMLAB, public grant
ANR-11-LABX-0037-01 constituting a part of the “Investissements
d'Avenir” program (reference : ANR-11-IDEX-0004-02; YK).






§ FITTING REPLICATION TIMING DATA FROM EXPERIMENTS USING THE
  MODEL


This section describes our fitting procedure based on the model. The
fitted parameters were used in simulations of genome replication
kinetics can giving the distribution of S-phase duration and of
replication time of one chromosome (Fig. <ref>).


We used flow cytometry (FACS) data to re-normalize replication timing
as follows.  If the base line value of average DNA copy-number a is
remarkably larger than 1, and/or its plateau value b is remarkably smaller
than 2, we use the formula
y=a+(b-a)(t-T_0))^r/(t-T_0)^r+(t_c-T_0)^rθ(t-T_0) to
fit the FACS data and normalize replication timing data by ϕ_
  norm(x,t)=1+ϕ(x,t)-a/b-a, where ϕ is the
replication probability function <cit.>.

We used fixed origin locations from the literature and optimized the
fit for the parameters γ, T_0, v and λ_i
iteratively. The objective function was defined as the L2 distance
(the average of squared differences) of the experimental and
theoretical replication probability timing profile (Fig. <ref>),
i.e., as √(∑_i∑_j(ϕ_model(x_i,t_j)-ϕ_exp.(x_i,t_j))^2/(N_xN_t)),
where N_x and N_t are the numbers of the measured loci and time
points respectively.

Initialization of the parameters for the fits was performed as
follows. Firing rate exponent γ and fork velocity v were
initialized at arbitrary values (typically γ at 0, v at 2
kb/min). The start of S phase T_0 was initially set when genome copy
number from the normalized FACs data (from the interval [a,b] to
[1,2]) is first larger than a fixed threshold (e.g. 1.05) and each origin
strength λ_i starts from the value fitted with the time-course
data at this origin.

Fitting was performed with following iterative rule. 1) for a
parameter x, assume it has a step length Δ_x, and a memorized
step length Δ_x^'=2Δ_x, 2) set r=Δ_x/Δ_x^'
and Δ_x^'=Δ_x, if x+Δ_x gives a better fit than
x, let x=x+Δ_x, otherwise (i) if |r|=1, we update
Δ_x →Δ_x/2 (ii) if |r|=0.5, set Δ_x
→ -Δ_x; 3)
repeat 2) until the termination condition is satisfied.  λ_1,
λ_2, ..., λ_n for each chromosome are updated
iteratively given γ, v and T_0 and in each iteration, one
λ_i is chosen randomly to be updated. T_0 is updated
iteratively given γ and v.  v is updated iteratively given
γ. For γ, we tested some discrete values between 0 and 3.
Supplementary Fig. <ref>a,b indicate the best fit value of γ
for S.cerevisiae and L.kluiveri,
and Supplementary Fig. <ref>c shows one example of the best fit.




§ ROLE OF CHROMOSOME BOUNDARIES IN REPLICATION TIMING



In some simulations, we used circularized chromosomes for easier
comparison with the analytical estimates, but relative to a circular
chromosome, a linear chromosome has lower symmetry because of the
boundary at both ends.  To verify that this assumption does not
qualitatively affect the results, we circularized the empirical
S.cerevisiae chromosomes by linking their ends respectively,
and simulated their replication kinetics with the estimated
parameters. The results (Fig. <ref>) show that the
circularized chromosomes always replicate faster than the linear
chromosomes, but their durations do not differ much (the average
deviation is in all cases less than 15%).










§ DETERMINATION OF THE PARAMETERS Α, Β AND T_0
  IN THE FORMULA FOR THE DISTRIBUTION OF T_I





Eq. 3 in the main text, describing the replication timing of one
inter-origin region contains the parameters α, β and
t_0, which need to be related to the biologically measurable
parameters (inter-origin distance and origin rates).  To estimate such
parameters for the distribution of T_i we used two methods. The
first is a fit of all the T_i data taken from the simulation of the
given chromosome, and the second is to fit the specific
T_i data (replication times of the central inter-origin region) extracted
 from simulation of a linear chromosomal fragment where inter-origin distances and origin
strengths are sampled from known distributions (different samples for different runs of the simulation).
In this second
method, each run of the simulation is carried out considering
inter-origin distances and origin strengths with the same averages as
the original chromosome.

Both methods give the same distribution for T_i, which agrees
very well with Eq. 3 of the main text (See Fig. <ref>).


We mainly used the second method since it does not depend on origin
configuration of the original chromosome. The detailed procedure is
the following. First, we defined a characteristic distance
d_c=(γ+1/<λ>log(1/1-x))^1/1+γv,
where x<1 (e.g. 0.99) and assume n_c=min(⌊
d_c/<d>⌋+1,⌊ n/2 ⌋)+1. Then we produced
a linear chromosomal fragment with 2n_c origins, in which two origins are
always located at the ends. Next, we simulated many realizations for
the replication of this chromosome. In each simulation run, we sampled
inter-origin distance d_i, origin strength λ_j and origin
firing time t_f^(j) from
Γ(<d>^2/σ^2(d),<d>/σ^2(d)),
Γ(<λ>^2/σ^2(λ),<λ>/σ^2(λ))
and f(t)=λ_i t^γθ(t)
  exp(-λ_it^γ+1/γ+1) respectively, where
i∈{1,2,...,2n_c-1} and j∈{1,2,...,2n_c}. The statistics
over different realizations gives the distribution of the replication
time of the central inter-origin region (T_n_c), which was fitted
with Eq. <ref> to obtain α, β and t_0.






§ ANALYTICAL DERIVATION OF AN APPROXIMATE DISTRIBUTION OF S-PHASE
  DURATION T_S BASED ON EXTREME VALUE THEORY.




This section gives further details on the analytical calculation for
the extreme-value estimate of the distribution of S-phase duration.
We assume that replication timing of one inter-origin region T_i
obeys the stretched exponential distribution

    F(t)=P(T_i<t)=1-e^-α (t-t_0)^β ,

where t⩾ t_0 and α>0. The parameters α,
β and t_0 were obtained as described in the previous
section. We define M_n=max(T_1, T_2,...,T_n). By
taking a_n=1/(α^1/ββ(log n)^1-1/β) and
b_n=(log n/α)^1/β+t_0, and applying the
Fisher-Tippett-Gnedenko theorem, we can prove that

    lim_n→∞P((M_n-b_n)/a_n≤
    t)=exp(-exp(-t))≜ G(t)  ,

where G(t) is the standard Gumbel distribution.

When n is sufficiently large, we can make the approximation
P((M_n-b_n)/a_n≤ t)≈ G(t).  If we define t̃=a_n
t+b_n, we have P(M_n≤t̃)≈ G((t̃-b_n)/a_n).

Finally, we can represent the distribution of T_S (=M_n)
approximately as

    P(T_S≤ t)  ≈exp(-exp(-t-b_n/a_n)) 
       = exp{-exp[βlog
      n(1-(α/log
      n)^1/β(t-t_0))]}








 Here n is the origin number, and α, β and t_0
  are connected to the model parameters describing replication
  kinetics, v, γ, inter-origin distances (d_1, d_2, ...,
  d_n) and origin strengths (λ_1,λ_2,...,λ_n).


We now discuss how α, β and t_0 can be expressed
  as functions of simplified parameters by numerically solving some
  approximate equations.  We consider a “characteristic”
  inter-origin region with the distance <d> and origin
  strength <λ>, and we assume that the replication
  of the inter-origin region is mainly carried out by the forks
  originated from the two nearest origins, both of which are typically
  activated, Thus we have

    T_i≈<d>/2v+(t_f^l+t_f^r)/2,

where t_f^l and t_f^r are the firing time of the left origin and
the right origin respectively. Since t_0 is the minimal replication time
of inter-origin region and the firing time has zero as a lower bound,
one has

    t_0=min(T_i)=<d>/2v.

From equation <ref>, we can further obtain

    <T_i> ≈<d>/2v+<t_f>

and

    σ(T_i)≈σ(t_f)

In addition, we have

    <T_i>=α^-1/βΓ(1/β+1)+t_0,


    σ(T_i) =
      α^-2/β[Γ(2/β+1) -
      Γ^2(1/β+1)],


    <t_f> =
      (γ+1/<λ>)^1/γ+1Γ(γ+2/γ+1),

and

    σ(t_f)=(γ+1/<λ>)^1/γ+1√(Γ(γ+3/γ+1)-Γ^2(γ+2/γ+1))

Based on equations <ref>-<ref>, α and β can be
numerically solved as functions of v, γ, <d> and
<λ>.  Our simulations in the EVD regime, and using
empirically realistic values of the parameters are in line with
equations <ref>-<ref>.





*













§ SUPPLEMENTARY FIGURES AND TABLES





































































































































































