BACKGROUND
the coordination of dna replication with cell division in e. coli is a classic problem of bacterial physiology  <cit> .
it is connected with the control of the bacterial dna replication and cell division cycle as a function of the growth rate, and it is an essential component for evolutionary adaptation to fast-growing conditions  <cit> .
it is also a classic problem for biological modeling  <cit> .
the main outstanding questions have to do with the characterization of the network of regulatory interactions by which cells determine the timing of initiation and limit it to once per cell cycle.
the theoretical foundations for understanding chromosome replication initiation in e. coli were set by cooper and helmstetter  <cit> , by showing that the time taken for a single chromosome to be replicated  and the time period between completion of chromosome replication and the following cell division  were approximately constant for a cell doubling time of less than one hour  <cit> .
the same work also introduced the idea of overlapping rounds of chromosome replication, where a round of replication can be initiated while an existing round of replication is still proceeding .
this mechanism allows e. coli to grow with a doubling time faster than the time required to copy its genome.
the question then arises of how the cell determines when to initiate dna replication and how this is coupled to the growth rate.
in  <dig>  donachie calculated that the correct timing would be guaranteed by a constant ratio of the cell size at the moment of initiation  and the number of oric in the cell  <cit> .
direct measurements of this ratio or of the initiation mass from cell population are difficult.
thus, whether this ratio is effectively a constant in cells with doubling times lower than one hour is to some extent an open question  <cit> .
it was then proposed and debated that the amount of an initiating factor accumulating with the cell's mass could reach a threshold value resulting in activation of the origin.
the dnaa protein has been shown to possess the basic characteristics necessary to act as such an initiator  <cit> .
several monomers of dnaa bind cooperatively to oric and induce dna melting required for assembly of the replication forks  <cit> .
its level of expression increases with the growth rate  <cit>  to result approximately in a constant amount of total dnaa per cell  <cit> .
dnaa overexpression results in earlier initiation, and when it is depleted it results in delayed initiation  <cit> .
in addition, dnaa exists under two forms, atp or adp bound.
the first is required for activation of the origin, thus it is usually called the active form  <cit> .
after the replication of oric, dnaa-atp becomes converted to dnaa-adp in a process known as 'rida' .
rida is mediated by the hda protein and the beta clamp subunit of the replisome and requires active replication forks  <cit> .
the hydrolysis of the atp in a dna-replication dependent manner decreases the activity of the protein after initiation has taken place, thus reducing the probability that a new initiation event will occur within the same cell cycle  <cit> .
at the same time the synthesis of new dna creates new dnaa binding sites that can titrate dnaa from the origin  <cit> .
other processes can contribute to prevent reinitiation within the same cell cycle, such as the binding of the seqa protein to the newly replicated, hemimethylated dna  <cit> .
it is believed that an "eclipse" period where reinitiation is not possible allows a buffer time for the other processes such as rida to take effect and thus for the levels of dnaa-atp to decrease below the critical level for initiation.
several gatc sites are also found at the dnaa promoter but their effect on the timing of initiation remains to be established  <cit> .
finally, a set of proteins have been shown to either inhibit or enhance the activity of dnaa at the origin.
these are for the most part abundant nucleoid proteins such as fis, hu and ihf, that may play a regulatory role as a function of changes in the growth phase  <cit> .
dnaa-atp binding to the origin must determine the timing of initiation for a range of growth rates and thus in the presence of increasing genome amounts .
thus, the amount of dnaa-atp per cell needs to increase with the decrease in doubling time.
the dnaa gene is found next to the origin on the chromosome, resulting in the gene copy-number increasing with the number of origins.
in addition, the expression of the dnaa gene is growth rate-dependent  <cit> .
the dnaa promoter region contains multiple binding sites for dnaa with differential affinity and specificity for the atp- and adp-bound forms of the protein and has been shown to be autorepressed by dnaa-atp but not dnaa-adp  <cit> .
consistent with this negative autoregulation, the artificial addition of dnaa boxes in the cell results in an increase in gene expression  <cit>  and inhibition of the rida process or the presence of a mutant form of dnaa insensitive to rida  results in a decrease in the level of dnaa protein in the cell  <cit> .
finally, dnaa is a transcription factor for a set of genes involved in regulation of dna replication  <cit>  and it could thus act as a reporter of the dna replication state of the cell in order to maintain the correct stoichiometry of the dna replication regulatory factors at varying growth rates and in response to perturbation to the movement of the replication forks  <cit> .
it has previously been proposed that the presence of both autoregulation and rida contributes to increased robustness of the initiation regulatory network upon perturbations  <cit> .
in this work, we aim to determine the relative roles of of these two regulatory processes in the control of the timing of initiation with changing growth rate.
we begin from the elements provided by the cooper and helmstetter model in order to estimate the initiation time at different growth rates.
the two main assumptions are that initiation of dna replication is determined by a critical amount of dnaa-atp per non-specific site on the genome and that this threshold value remains constant as a function of growth rate.
on the other hand, the cellular and metabolic parameters can change with growth rate and have an impact on the dnaa circuit.
in order to understand this, we use information from systematic studies of cellular changes with growth rate  <cit> .
finally, the volume of the cell is assumed to be a less relevant background as a reservoir of dnaa-atp than the number of non-specific binding sites on the dna  <cit> .
the resulting equations describe, via a continuous change in parameter values with growth rate, the oscillations in dnaa-atp per non-specific site and the attainment of a constant threshold as a function of growth rate.
this shows that the circuit performing the timing of replication initiation must encode subtle information on the bacterial physiological state through the growth rate dependence of the parameters.
this analysis also allows us to define a few scenarios consistent with the available experimental knowledge and to make testable predictions on the relative roles of dnaa autorepression and of the rida process at different growth rates.
we use this model to elucidate the reciprocal roles of the known factors affecting dnaa activity in e. coli, namely that dnaa expression is dependent on transcriptional autoregulation, and that its atp-ase activity is coupled with the activity of the advancing replication forks .
the results show that a working system can still be produced in the absence of rida or dnaa autoregulation.
moreover, both rida and autorepression contribute to a decrease in the amplitude of the cell-cycle oscillations in dnaa-atp.
rida has a larger effect at the faster growth rates while negative regulation has a larger effect at slower growth.
methods
assumptions of the model
the model consists of a set of ordinary differential equations  describing dnaa-atp production by the expression of the dnaa gene.
it is built on two basic assumptions.
the first is based on the evidence that a specific number of dnaa-atp molecules need to bind on oric in order to create a replication bubble.
following the standard thermodynamic model of protein-dna binding  <cit> , the probability of this event is dependent on the number of dnaa-atp molecules that are bound to the non-specific sites along the chromosome.
these are low-affinity sites compared to titration sites, but the affinity is high enough so that the protein spends most of its time bound to the genome.
these sequence-independent interactions are typical of dna-binding proteins.
as a consequence, the simplifying assumption is usually made  <cit>  that the key molecular players  are bound to the dna either specifically or non-specifically.
simply stated, this is just an implementation of the known fact that dna-binding proteins, besides binding tightly to their target sequences, are generally "sticky" for dna, in a sequence-independent manner.
this implies that the timing of dna replication initiation in the cell is determined by the ratio dnaa-atp to non-specific binding sites on dna  in the cell and effectively results in the computation of the amount of dna per cell).
note that in this case the volume of the cell is not taken into consideration since it is not the change in concentration of dnaa or of dna that determines the initiation time.
the same will apply to the binding of dnaa and rna polymerase to the dnaa promoter .
the second assumption is that at the time of initiation this ratio will be the same, independently of the growth rate and thus the number of origins.
we ask how the parameters of this model must vary in order for this assumption to hold in the range of doubling times between  <dig> and  <dig> minutes.
in the absence of autoregulation, the only factor that contributes to a decrease in the ratio of dnaa-atp to non-specific binding sites is the increase in dna after dna replication has begun.
the complete model also includes the autoregulation of dnaa expression by dnaa binding to its own promoter  <cit>  and dnaa-atp transformation into dnaa-adp through the rida process  <cit> .
it is assumed that this ode description is applicable to a single average-cell on time scales shorter than the length of the cell cycle.
this hypothesis could be challenged for the shortest observable doubling times, but the formulation of the model is dictated by maximizing simplicity.
the values of the parameters  are all taken or estimated from the available experimental measurements.
they are shown in table  <dig>  together with the sources.
the parameters are fixed with the values in the table for τ =  <dig> mins and then are able to vary for the other values of τ in order to fix the ratio of dnaa-atp:chromosome length at the moment of initiation .
†: chiara sag-gioro, anne olliver, bianca sclavi: multiple levels of regulation in the growth rate dependence of dnaa expression, submitted.
formulation of the model
timing of replication
we take into consideration the situation where the cell cycle repeats itself identically i.e.
balanced, exponential growth.
following cooper and helmstetter, at a time c + d after the initiation time, the cell divides, i.e.
that time must be an integer multiple of the doubling time.
thus, if τ is the doubling time of the cell and x is the initiation time, then we must have
x+c+d=τ,
where n is the integer number of times that τ divides c + d. n can be viewed as the number of overlapping rounds of replication, and 2n is the number of origins.
thus, this equation reflects the phenomenon of overlapping replication rounds.
figure  <dig> shows how x varies with the doubling time  of the cell.
defining y as the time at which the chromosome completes replication, we have
y=τ-d.
promoter term
the activity of the dnaa promoter is the source term for dnaa-atp.
we describe it by the standard thermodynamic model first used by shea and ackers  <cit> .
we denote the promoter term  as q, the number of non-specific binding sites on the chromosome as nns and the number of rna polymerase  molecules as p. furthermore, we denote the number of dnaa-atp molecules as a-.
we then use the assumption that the number of non-specific binding sites is proportional to the length of dna in the cell , i.e.
Λ = κnns.
thus the expression for the rate of transcription at the dnaa promoter can be written as
q=kaΘ1+c1Λp+c2a-p
where ka is the basal rate of transcription of the dnaa promoter, Θ is the number of dnaa promoters  in a given cell at a given time, and the remaining factor is the probability of rnap binding to the promoter.
the parameters c <dig> and c <dig> depend on the binding energies Δεpd and Δεad of rnap and a- respectively to their promoter binding sites.
the binding energies are determined from the ratio of specific vs non-specific binding affinities.
c1=eΔεpd∕kbtκc2=e∕kbt.
where the exponential terms are boltzmann weights.
c <dig> =  <dig> if the promoter is not autorepressed.
a version of the promoter where dnaa binding to its sites is cooperative is described in additional file  <dig>
rida term
this term reflects the number of dnaa-atp molecules that are converted to dnaa-adp molecules per unit time by the rida process.
as discussed in the introduction, rida is a process that takes place at the replication forks during dna synthesis.
we assume that the rate of conversion kr takes the same value at each replication fork.
the number of pairs of replication forks at a given time, f, depends on which of x and y is larger.
for x < y :
f=2n-12⋅2n-12if0<t<xifx<t<yify<t<τ
and for x > y :
f=2n-12n-22if0<t<yify<t<xifx<t<τ
and so the conversion from dnaa-atp to dnaa-adp takes place at a rate
krf.
this leads to the following differential equations for dnaa-atp  and dnaa-adp
∂a-∂t=kaΘ1+c1Λp+c2a-p-krf
∂a+∂t=krf.
term for the growth of the chromosome
the growth of the chromosome is controlled by the replication forks.
defining the rate of dna synthesis of each pair of replication forks as kΛ, we can write
∂Λ∂t=kΛf.
assuming that kΛ is constant, and normalizing so that Λ =  <dig> is the length of one full chromosome, we have kΛ = 1/c.
main equation
∂r∂t=1ΛΘka1+c1Λp+c2Λpr-f.
this equation describes the dynamics of the variable r=a-Λ which we suggest is a suitable candidate for the initiation potential since a specific number of dnaa-atp molecules is needed to be available to bind to the origin in order to induce dna melting, as described above.
note that usually such dynamic equations are written in terms of volume, thinking of averages over cell populations on time-scales longer than a cell cycle.
we assume that the  background of genomic binding sites is the relevant variable in a shea-ackers type model, extending to a single cell cycle the approach normally used for longer time scales  <cit> .
we also assume that the volume can be treated as a weak perturbation, which we neglect here.
in other words, the various molecules of interest  are partitioned between the specific and non-specific binding sites on the chromosome.
furthermore, the most important factor that determines the probability of binding to a given promoter is the absolute number of protein molecules relative to the absolute number of these binding sites, rather than the amount of protein per cell volume, which in comparison does not change significantly, and it can thus be neglected to a good approximation.
thus, this assumption means that one need not track the volume of the cell, only the number of non-specific binding sites in the cell at a given time.
this idea is discussed further in additional file  <dig>  section  <dig>  here we assume that the initiation potential, r, always reaches the same value at t = x independently of the growth rate and we ask how the parameters of this model must vary in order for this assumption to hold in the range of doubling times between  <dig> and  <dig> minutes.
in the following, we will first establish that such a constant threshold cannot be obtained by a model with fixed parameters and then study the possible scenarios where different subsets of parameters are allowed to vary.
main assumptions
the model relies on the following further assumptions  <cit> :  all newly synthesised dnaa is immediately bound to atp, due to the relative abundance of atp in the cell compared to adp and the high affinity of dnaa for atp.
dnaa-adp is only created by conversion from dnaa-atp by the rida process, when it is present.
the probability of dnaa being bound to its sites at the origin or on the promoter, in the case of the presence of autorepression, is given by its thermodynamic equilibrium value.
the same assumption holds for the binding of rna polymerase to a particular promoter.
this means that we assume that the rate for transcription initiation is much slower than the rates for rnap binding and unbinding from the promoter.
the rate of dnaa gene expression is proportional to the equilibrium probability that rnap is bound to the dnaa promoter.
we do not consider translation directly and thus there is no time delay from transcription to protein production since the addition of this feature did not affect the result of the model .
the number of non-specific binding sites on the dna in each cell for both rnap and dnaa is proportional to the total length of dna in the cell  <cit> .
the number of rnap molecules in the cell, p, grows exponentially from cell birth to cell division, corresponding to the hypothesis of constant concentration and exponentially growing cell size  <cit> .
a linear growth can also be used, however the dynamics of the model do not differ significantly between these two cases.
numerical integration
the non-linearity of the main equation  necessitates the use of a numerical method of integration.
we used a custom c++ implementation of the fourth-order runge-kutta method.
the equation was integrated for values of the cell doubling time, τ, in the range  <dig> mins ≤ τ ≤  <dig> mins.
in order to test for the constant threshold condition, a transformation was performed by integrating the equation for τ =  <dig> mins and using the value of r at initiation t = x as the imposed threshold for the other values of the doubling time.
thus it was important to estimate, to a good degree of accuracy, values for the parameters at a doubling time of  <dig> mins.
these values appear in table  <dig>  the parameter values are either taken to be constant  or are allowed to change and obtained as a consequence of the transformation .
in all the scenarios, ka varies with growth rate.
in scenario 1a, p changes as a function of growth rate with the values obtained from ref.
<cit> .
the dependence of the other parameters on growth rate is scenario specific .
RESULTS
a fixed set of parameters gives a varying initiation threshold with increasing growth rate
we first describe the behaviour of the model with a fixed parameter set.
the ratio r=a-Λ increases from the time of birth of the new cell.
this can be interpreted as the accumulation of the 'initiation potential'.
at initiation , r peaks  and then falls again due to the increase in the number of non-specific binding sites.
when including the rida process, the total rida rate also increases following initiation, due to the higher number of active replication forks, contributing to the decrease in the initiation potential.
however, the value of r at t = x varies for the different values of the doubling time t .
thus, it appears that this model, with fixed parameters, cannot give a constant threshold that is reached at initiation in the range of growth rates considered here.
this fact naturally leads to us consider a model in which the parameters are able to vary with growth rate.
biologically, this is a natural requirement, as one may well expect from previous observations of the change in cellular components as a function of growth rate  <cit> .
a constant threshold condition implies alternative scenarios of growth-rate dependency in the circuit architecture
the condition of a constant dnaa-atp/dna threshold at the time of initiation can be imposed by performing a mathematical transformation on the model and verifying the implications of this for the values of the parameters.
the mathematical details of this transformation can be found in additional file  <dig>  section  <dig>  the transformation yields a fixed threshold in r that is reached at initiation.
this can be seen upon comparing the plots in figure 3a and 3b.
in the latter, the value of r is now the same at initiation for every cell doubling time.
we estimate all initial values of the parameters from the literature  and then we allow some of the parameters to change during the transformation.
thus, the transformation procedure imposes a decision on which parameters to fix and which to allow to change.
this determines different scenarios, as illustrated in figures  <dig> and  <dig> and also summarised in table  <dig> .
these plots explicitly show the required parameter changes at varying growth rates.
<dig>  in the first scenario, the rida rate , kr, is chosen to remain independent of growth rate, and the other parameters are allowed to change.
the transformation fixes the scaling of ka with growth rate, and provides an equation that c <dig> and p must satisfy.
since the constraints reduce to one equation with two unknowns, there is a set of possible solutions for c <dig> and p, and the scenario is underconstrained .
one possibility is to fix one of the parameters to a particular trend, in turn determining the second parameter, determining different subscenarios .
in the first of these sub-scenarios, 1a, the variation in the number of non-specifically bound rnap in the cell, p, with growth rate is fixed a priori by imposing the trend of a previous study  <cit> , which partitions rnap into different classes, the rnap bound on promoters and non-specifically bound, and determines the dependence of this partition with growth rate, .
imposing this trend determines the values for the binding affinity of rnap to the promoter, which must vary with growth rate .
in sub-scenario 1b, the binding affinity of rnap to the dnaa promoter is chosen to remain independent of growth rate.
fixing this parameter constrains the remaining equation, thereby imposing a particular dependence on growth rate of the number of non-specifically bound rnap, p, which turns out to be compatible with ref.
<cit>  .
<dig>  in scenario  <dig>  the binding affinity of both rnap to the dnaa promoter and of dnaa-atp to its repression sites are chosen to remain independent of growth rate , while the basal rate of transcription from the dnaa promoter, ka, the levels of free rnap, p, and the rida rate , kr are all permitted to vary from their original values .
this scenario is also feasible in absence of autorepression.
in brief, two possible scenarios can result in a constant initiation threshold for the model.
in the first, the binding affinity of dnaa-atp to its repressor sites decreases with increasing cell doubling time and in the second the rida rate increases with cell doubling time.
in both scenarios ka  must decrease with increasing cell doubling time.
it is important to note that ka, the basal dnaa transcription rate, needs to vary in all scenarios with growth rate.
this term  is independent from rnap availability  and binding .
it describes how quickly rnap moves through a gene when transcribing.
the variation of this characteristic time with growth rate can be associated with variations in dna supercoling .
in absence of rida  or autorepression , the transformation can still be performed.
however, it implies that the ratio c1/p should remain constant with growth rate.
since p varies  <cit> , this means that c <dig>  would have to compensate exactly for the changes in p in order to keep a constant threshold.
we have considered these further scenarios to be less probable, because they may result in a less robust control due to an unlikely fine-tuning of two parameters.
the resulting scenarios are compatible with available knowledge on rnap availability and total dnaa expression
given these scenarios, we have asked whether the predicted parameter variation with growth rate and the observables quantities produced by the model were compatible with the measurements and observations available in the literature.
starting from the dependency of available rnap with growth rate, this is predicted and matched with available experimental data in ref.
<cit> .
scenario 1a assumes this dependency, and therefore automatically accounts for this observation.
figure 4a demonstrates how also scenarios 1b and  <dig> are broadly compatible with the results of klumpp and hwa  <cit>  since the average levels of non-specifically bound rna polymerase decrease with increasing doubling time in all cases.
we now turn to the changes in measured expression of dnaa  with growth rate.
this can be measured by a reporter gene technique.
figure 4b shows how the model predicts that the average expression of the dnaa gene should change with growth rate in all scenarios.
this appears to be independent of the scenario chosen and is compatible with previous findings  <cit> .
we have also performed our own measurements, using a gfp reporter of the dnaa promoter encoded on a plasmid , see also additional file  <dig>  section 7), confirming this agreement .
finally, at fixed growth rate, in order to determine whether this model reflects the main features of the regulatory network in the cell, we reproduced some of the experimental perturbations described in the literature.
one of these experiments changed the rate of rida by changing the level of expression of the gene encoding for the hda protein  <cit> , while others controlled the expression of dnaa independently of the dnaa promoter.
we found that a constant threshold can be obtained upon a  <dig> fold change in the rida rate .
in order to determine how rida rate and autorepression strength influence the activity of dnaa-atp within a cell cycle in the model, we have monitored the amplitude of the oscillations of the ratio dnaa-atp:total genome length.
the results show that increased autorepression contributes to a smaller amplitude of dnaa-atp within the cell cycle .
a decrease in rida has the opposite effect, and is compensated for by an increase in autorepression.
in either case the inclusion of these additional control factors appears to result in smaller amplitude in the oscillation of dnaa-atp and a smaller variation in the average amount of dnaa-atp per cell as the growth rate is varied.
this may be advantageous for the use of dnaa as a transcription factor whose activity is responsive to changes in the replication status of the cell, via the rida and titration effects.
finally, once we obtained a set of parameters that satisfies the constant threshold constraint, we modified the rida rate while leaving the other parameters unchanged , attempting to reproduce the effect of under or over expressing the hda protein, as in the experiments by riber et al.
<cit> .
in our model, this results in a change of the threshold value as a function of growth rate, and more significantly at slow growth.
at faster rida rate, the threshold value is higher at slow growth, while the opposite is observed when rida rate is decreased.
these compatibility tests give positive results, but do not allow us to distinguish between the two scenarios.
we have explored the literature for tests of dependency of the rida rate with growth rate, and have found no evidence of this, which lead us to consider scenario  <dig>  where the rida rate is constant as the main one.
model variants
in order to gain confidence that the conclusion  is not a consequence of the restricted set of biological ingredients included in this model, we considered some additional model variants, including some of the known factors that may influence the timing of replication initiation.
<dig>  delay in the synthesis of dnaa-atp
we introduced a delay, representing the time necessary to obtain an active dnaa molecule from the binding of the rna polymerase to the dnaa promoter to the end of translation.
this delay, however, does not produce a significant effect in imposing a dnaa-atp threshold at initiation, suggesting that translation delay might not have a predominant role in controlling the timing of initiation .
<dig>  cooperativity of autorepression
cooperativity of autorepression affects the growth rate dependence of gene expression  <cit> .
additional file  <dig>  figure a6c shows by simulation that the presence of cooperativity alone cannot explain a constant threshold.
thus it is still necessary to impose a constant threshold on the model using a transformation such as that described in additional file  <dig>  section  <dig>  however, in general one cannot make a translation and scaling and keep the promoter in the same form , and the transformation itself poses complex constraints on the possible regulations.
thus, we decided to concentrate on the simpler non-cooperative model.
biologically, of the two high-affinity dnaa sites found at the dnaa promoter, one matches exactly the consensus sequence, and has a higher affinity than the other  <cit> , suggesting that at lower dnaa concentrations only one monomer could be bound.
we have verified that the transformation is possible with a promoter that would fit this profile, i.e.
of the form
q=kaΘ1+c1Λp1+rk1+r2ωk1k21+rk <dig>
where the parameter ω represents the cooperativity, and where k <dig> and k <dig> are the binding affinities of the two dnaa binding sites, multiplied by the proportionality constant between Λ and nns.
this description introduces two new parameters - the cooperativity and the binding affinity of the second binding site.
since these parameters appear in the equations describing the transformation that keeps the constant threshold, the scenarios become more underconstrained in a non-essential way.
thus, while a more realistic promoter model could be useful for future descriptions, we decided not to pursue it here.
<dig>  the data locus and specific dnaa binding sites
we considered the effect that the presence of the data locus has on the model by introducing a site on the chromosome that binds up to  <dig> dnaa-atp molecules immediately after initiation has taken place .
we included one data site for each origin in the cell.
the results indicate that the data locus might indeed prevent further initiations in a given cell cycle, since it titrates large numbers of dnaa-atp molecules, effectively preventing them from binding at the origin .
however, this variant of the model fails to achieve a constant dnaa-atp threshold at initiation at different growth rates if all parameters are kept constant.
one can speculate that the large increase in data sites at faster growth rates would require a proportional increase in the rate of dnaa synthesis that cannot be solely provided by the increase in gene copy number.
it must be also noted that the data locus appears to be unnecessary to prevent reinitiation and limit initiation of to once per cell cycle  <cit> .
additionally, we considered the effect of the reported ≈  <dig> binding sites distributed around the chromosome  <cit> .
these high-affinity sites would sequester dnaa in a similar way as the data locus, but proportionally to the genome amount, and thus effectively decrease its concentration.
this is again insufficient to guarantee a constant initiation threshold and is equivalent to rescaling the rida rate .
<dig>  the eclipse
no constraint for the eclipse period  <cit>  was included in this model.
this proved not to be necessary, probably because there is no delay between the attainment of the threshold and the initiation of dna replication, thus avoiding the possibility that there is an 'overshoot' of initiation potential before both rida and the increase of non-specific dna can begin.
on the other hand, the gene copy number immediately doubles upon initiation of dna replication which results in a sudden increase in the number of dnaa-atp synthesized per unit time.
however, this does not result in an increase in the initiation potential due to the corresponding increase in the number of replication forks and thus on the rida rate and the number of non-specific titrating sites.
<dig>  dnaa-atp recycling regions
genomic recycling regions catalyzing the reconversion of dnaa-adp into dnaa-atp  <cit>  have been recently discovered.
the quantitative contribution of this process is not clear, but within our framework it makes sense to ask how this would affect the initiation threshold, as this process is, roughly, a correction to rida.
precisely, assuming dnaa-adp is not limiting, this would change the model by the substitution of kr with kr-ρΛ∕f where ρ is a fixed recycling rate.
we verified that in the model this recycling term cannot by itself impose a constant threshold, while it can contribute to correcting quantitatively the effective rida rate .
since none of these model variants qualitatively changes the behaviour of the model with respect to attaining a constant initiation threshold, they were not included in the minimal model formulation, in order to avoid confusion and proliferation of parameters.
however, as shown by the variants explored above, the qualitative behaviour that the parameters of the model must vary with growth rate does not hold strictly for the minimal model only, but might be more general.
discussion
standard models of bacterial regulatory circuits were adapted to situations where the growth rate is fixed  <cit> .
the notion that these quantitative descriptions must account for bacterial physiology through the growth-rate dependent basic partitioning of the cell physico-chemical components is now entering the field of systems biology through a combination of new work  <cit>  and reconsideration of the classics  <cit> .
the dependency of the basic parameters on growth rate can produce notable effects on a genetic circuit, and complicates the standard descriptions  <cit> .
in our case, the task is more difficult, as the circuit under examination is active in determining some features of the bacterial physiology and not only affected by them.
furthermore, on the technical level, one must produce a time-dependent description the expression of dnaa over cell-cycles of a range of durations.
perhaps also for this reason, despite the fact that the regulation of dna replication has been a subject of intense study for over  <dig> years  <cit> , many questions remain open.
given these obstacles, we have shown that, under a series of simplifying hypotheses, a consistent mean-field description for the dnaa/replication initation circuit is possible with varying growth rate.
our description includes the processes that are believed to be most important for initiation of replication  <cit> .
in these respects, it is broadly compatible with previous modelling approaches  <cit> .
its originality lies in the minimality and in the attention given to growth-rate dependency.
we focused on the minimal ingredients necessary in order for the basic tenet that the ratio dnaa-atp/dna attains a constant threshold at initiation to hold  <cit> .
the validity of this tenet is confirmed by the recent observations that initiation time is not affected by adding an extra origin on the chromosome  <cit>  and on the compensatory mutations emerging in hda mutants  <cit> .
we have defined the dna replication initiation potential, determining the  timing of dna replication, as the dnaa-atp to dna ratio, r. molecular titration has been shown to result in ultrasensitive "all or none" responses  <cit> , which further justifies using r as the threshold and could explain the synchrony of initiation in cells containing oric minichromosomes  <cit> .
we assume that its value at the time of initiation, r, is independent of the specific growth rate.
the amount of dnaa-atp at the time of initiation thus needs to increase as a function of growth rate in order for r to remain constant as a function of doubling time, and we found that consequently, some of the model's parameter values must be allowed to vary.
this assumption has not been verified directly.
on the other hand, we feel that our point of view would still be useful in case of a growth-rate dependent r, as it is unlikely that this dependence would automatically match the dependence of all the other parameters.
we have defined two main scenarios in which different subsets of the parameters are allowed to change.
in scenario  <dig> the rida rate , kr, is held constant as a function of growth rate, but the binding affinities of rnap and dnaa-atp to the dna need to vary with growth rate  or allowing it to be free ).
in scenario  <dig> the binding constants  are independent of growth rate but the rida rate, kr, must vary.
we have verified that both scenarios are consistent with the eperimentally tested predictions of rnap availability with growth rate  <cit>  and with previous measurements  <cit>  and our own experimental evaluation of total dnaa expression , and also with a number of "in silico mutations" inspired by the available literature  <cit> .
thus, the scenarios appear as possibilities that are testable, but for the moment remain open.
note that the property that the initiation threshold holds constant with respect to growth rate changes is not related to the specific set of parameters we used, or any set of parameters.
our analysis shows that in general, for any fixed parameter set at a given growth rate, a transformation is necessary in order to keep the threshold constant while moving to another growth rate.
in order to provide specific examples, we have produced plots in the style of those in figure  <dig>  with different curves corresponding to choosing different values of the initial input parameters.
these demonstrate that the qualitative behaviour of the transformation is independent of these parameters .
this exercise is also important to show that the parameter changes with growth rate are not numerically negligible for empirically plausible parameters, so that the question of keeping the initiation threshold constant is not purely academic.
it is then interesting to ask which of these scenarios is more reasonable considering the known biological processes.
we speculate that scenario  <dig> is less likely, since until now there is no evidence pointing to a possible change in the intrinsic rida rate as a function of growth rate.
the dnaa-related protein hda  mediates this process  <cit> .
experiments with mutants over-and underexpressing hda  <cit> , with corresponding increases and decreases in the rida rate, suggest a possible mechanism by which the kr term in the equations could vary by a growth rate-dependent expression of the hda gene.
there may also be other, as yet unknown, factors that affect the growth rate dependence of the rida rate.
alternatively, we can speculate that the decrease in the rate of rida with growth rate could be caused effectively by the action of the reverse process of dnaa-atp recycling by the recently discovered recycling regions  <cit> .
figure 5b shows that the rida rate should increase with cell cycle time and thus decrease with growth rate.
this growth rate increase causes overlapping replication rounds, and thus higher chromosome copy number.
since more recycling regions are present there is more recycling, i.e.
a decrease in the effective rida rate, compatibility with the requirement imposed by our results.
however, considering explicitly this model variant, we find that the balancing recycling cannot by itself impose a constant threshold.
conversely in scenario  <dig>  the rida rate per replication fork is constant, and one has to rationalize the variation of the binding affinities.
it seems possible that the binding affinities could change with growth rate through changes in supercoiling, in similar ways to those seen in figure  <dig> and additional file  <dig>  figure a <dig>  <cit> .
the levels of average negative supercoiling are known to increase as the growth rate increases  <cit> .
however, it makes sense to challenge the validity of the basic assumption that the ratio of dnaa-atp to dna at the time of initiation is constant.
this model assumes that the affinity for dnaa-atp binding to its own promoter can change with growth rate but its affinity for the origin does not.
the first assumption mainly allows the model to change the magnitude of negative autoregulation as a function of growth rate, and it may indeed be explained by the changes in global cellular parameters such as negative supercoiling.
we have considered how the activation threshold in the model  would be affected if the binding affinity for dnaa to the origin would vary in the same way as its value at the dnaa promoter, required by scenario  <dig>  for a set of realistic parameters.
we found that these changes in r are less than 10% over a wide range of growth rates, suggesting that this scenario might be robust.
indeed, the observed threshold is certainly approximately constant when compared to the untransformed case i.e.
the different values of the ratio at t = x shown in figure 3a.
more generally, the initiation of dna replication has been significantly simplified in this model; all it requires is a specific amount of available dnaa-atp molecules.
however we know that other factors, such as the binding of nucleoid proteins fis, ihf, h-ns and hu, may contribute to the formation of an open complex at the origin.
on the other hand, other recent results have shown that at slower growth  the cell contains a greater average amount of dnaa-atp per origin that results in initiation events that are independent of the novel synthesis of dnaa-atp  <cit> .
these results suggest that the regulation of the initiation process at the origin might indeed be dependent on the growth rate and that these changes still remain to be characterized quantitatively before they can be included in a theoretical model.
interestingly, the basal rate of transcription of the dnaa gene, ka must vary in both scenarios.
figure  <dig> shows that ka decreases as the cells grow more slowly.
this is what is expected from a promoter like the one of the dnaa gene that closely resembles ribosomal rna promoters.
this family of promoters have a gc-rich sequence at the transcription initiation site called a discriminator region.
this region renders the activity of the promoter sensitive to the degree of negative supercoiling, which activates transcription by enhancing dna melting, and leads to its inhibition by the accumulation of ppgpp at slower growth rates  <cit> .
CONCLUSIONS
all things considered, we can say that perhaps our main result is that the determination of the timing of initiation by dnaa, besides relying on the known "architecture" comprising autorepression, rida and a number of other "dedicated" processes, can be understood only in its complex interplay with bacterial physiology
nevertheless, it makes sense to ask whether this model allows us to elucidate some features of the reciprocal role of rida and dnaa autorepression, its two main ingredients.
biologically, rida renders the control of dnaa-atp dependent upon ongoing dna replication, and thus results in an increase in dnaa-atp when replication forks are blocked.
autorepression however probably plays a larger role in the absence of rida at slow growth, or in bacteria that do not have rida at all   <cit> .
the hda protein and thus the rida process seems to be quite specific to the fast-growing e. coli bacterium and its close relatives in the enterobacteriaceae , suggesting that in other bacterial species this level of regulation may not be required and is replaced instead by protein degradation, e.g.
in caulobacter  <cit> , or a high intrinsic atpase activity of the protein, as in mycobacterium tuberculosis  <cit> .
we have verified that the model can work in the absence of autorepression or rida, but the tuning of the parameters to achieve a constant threshold is more "difficult", in the sense that it requires more fine-tuning of the parameters, since the ratio c1/p should remain constant with growth rate.
meaning that it is possible that a smaller range of growth rates would be accessible in these conditions.
moreover, in the model, increasing autorepression or rida rate results in a smaller amplitude of the oscillations of the ratio dnaa-atp/dna during the cell cycle, and in a smaller variation in the average amount of dnaa-atp per cell as the growth rate is varied.
this may be advantageous for the use of dnaa as a transcription factor which has to sense perturbations in the replication status of the cell at all growth rates.
competing interests
the authors declare that they have no competing interests.
authors' contributions
mcl, bs, and bb designed research.
mg, bb, mcl, uf, and cs performed research.
mg, bs and mcl wrote the paper.
all authors read and approved the final manuscript.
supplementary material
additional file 1
additional text and figures.
single pdf file containing the additional text and figures mentioned in the main text.
click here for file
acknowledgements
we are grateful to matteo osella, rosalind allen, pietro cicuta, antonio celani, andrea sportiello, kunihiko kaneko and massimo vergassola for useful discussions and feedback.
the authors acknowledge support from the human frontier science program organization  and from epsrc.

 the coordination of dna replication with cell division in e. coli is a classic problem of bacterial physiology. it is connected with the control of the bacterial dna replication and cell division cycle as a function of the growth rate, and it is an essential component for evolutionary adaptation to fast-growing conditions. it is also a classic problem for biological modeling. the main outstanding questions have to do with the characterization of the network of regulatory interactions by which cells determine the timing of initiation and limit it to once per cell cycle. the theoretical foundations for understanding chromosome replication initiation in e. coli were set by cooper and helmstetter, by showing that the time taken for a single chromosome to be replicated and the time period between completion of chromosome replication and the following cell division were approximately constant for a cell doubling time of less than one hour. the same work also introduced the idea of overlapping rounds of chromosome replication, where a round of replication can be initiated while an existing round of replication is still proceeding. this mechanism allows e. coli to grow with a doubling time faster than the time required to copy its genome. the question then arises of how the cell determines when to initiate dna replication and how this is coupled to the growth rate. in donachie calculated that the correct timing would be guaranteed by a constant ratio of the cell size at the moment of initiation and the number of oric in the cell. direct measurements of this ratio or of the initiation mass from cell population are difficult. thus, whether this ratio is effectively a constant in cells with doubling times lower than one hour is to some extent an open question. it was then proposed and debated that the amount of an initiating factor accumulating with the cell's mass could reach a threshold value resulting in activation of the origin. the dnaa protein has been shown to possess the basic characteristics necessary to act as such an initiator. several monomers of dnaa bind cooperatively to oric and induce dna melting required for assembly of the replication forks. its level of expression increases with the growth rate to result approximately in a constant amount of total dnaa per cell. dnaa overexpression results in earlier initiation, and when it is depleted it results in delayed initiation. in addition, dnaa exists under two forms, atp or adp bound. the first is required for activation of the origin, thus it is usually called the active form. after the replication of oric, dnaa-atp becomes converted to dnaa-adp in a process known as'rida'. rida is mediated by the hda protein and the beta clamp subunit of the replisome and requires active replication forks. the hydrolysis of the atp in a dna-replication dependent manner decreases the activity of the protein after initiation has taken place, thus reducing the probability that a new initiation event will occur within the same cell cycle. at the same time the synthesis of new dna creates new dnaa binding sites that can titrate dnaa from the origin. other processes can contribute to prevent reinitiation within the same cell cycle, such as the binding of the seqa protein to the newly replicated, hemimethylated dna. it is believed that an "eclipse" period where reinitiation is not possible allows a buffer time for the other processes such as rida to take effect and thus for the levels of dnaa-atp to decrease below the critical level for initiation. several gatc sites are also found at the dnaa promoter but their effect on the timing of initiation remains to be established. finally, a set of proteins have been shown to either inhibit or enhance the activity of dnaa at the origin. these are for the most part abundant nucleoid proteins such as fis, hu and ihf, that may play a regulatory role as a function of changes in the growth phase. dnaa-atp binding to the origin must determine the timing of initiation for a range of growth rates and thus in the presence of increasing genome amounts. thus, the amount of dnaa-atp per cell needs to increase with the decrease in doubling time. the dnaa gene is found next to the origin on the chromosome, resulting in the gene copy-number increasing with the number of origins. in addition, the expression of the dnaa gene is growth rate-dependent. the dnaa promoter region contains multiple binding sites for dnaa with differential affinity and specificity for the atp- and adp-bound forms of the protein and has been shown to be autorepressed by dnaa-atp but not dnaa-adp. consistent with this negative autoregulation, the artificial addition of dnaa boxes in the cell results in an increase in gene expression and inhibition of the rida process or the presence of a mutant form of dnaa insensitive to rida results in a decrease in the level of dnaa protein in the cell. finally, dnaa is a transcription factor for a set of genes involved in regulation of dna replication and it could thus act as a reporter of the dna replication state of the cell in order to maintain the correct stoichiometry of the dna replication regulatory factors at varying growth rates and in response to perturbation to the movement of the replication forks. it has previously been proposed that the presence of both autoregulation and rida contributes to increased robustness of the initiation regulatory network upon perturbations. in this work, we aim to determine the relative roles of of these two regulatory processes in the control of the timing of initiation with changing growth rate. we begin from the elements provided by the cooper and helmstetter model in order to estimate the initiation time at different growth rates. the two main assumptions are that initiation of dna replication is determined by a critical amount of dnaa-atp per non-specific site on the genome and that this threshold value remains constant as a function of growth rate. on the other hand, the cellular and metabolic parameters can change with growth rate and have an impact on the dnaa circuit. in order to understand this, we use information from systematic studies of cellular changes with growth rate. finally, the volume of the cell is assumed to be a less relevant background as a reservoir of dnaa-atp than the number of non-specific binding sites on the dna. the resulting equations describe, via a continuous change in parameter values with growth rate, the oscillations in dnaa-atp per non-specific site and the attainment of a constant threshold as a function of growth rate. this shows that the circuit performing the timing of replication initiation must encode subtle information on the bacterial physiological state through the growth rate dependence of the parameters. this analysis also allows us to define a few scenarios consistent with the available experimental knowledge and to make testable predictions on the relative roles of dnaa autorepression and of the rida process at different growth rates. we use this model to elucidate the reciprocal roles of the known factors affecting dnaa activity in e. coli, namely that dnaa expression is dependent on transcriptional autoregulation, and that its atp-ase activity is coupled with the activity of the advancing replication forks. the results show that a working system can still be produced in the absence of rida or dnaa autoregulation. moreover, both rida and autorepression contribute to a decrease in the amplitude of the cell-cycle oscillations in dnaa-atp. rida has a larger effect at the faster growth rates while negative regulation has a larger effect at slower growth. methods assumptions of the model the model consists of a set of ordinary differential equations describing dnaa-atp production by the expression of the dnaa gene. it is built on two basic assumptions. the first is based on the evidence that a specific number of dnaa-atp molecules need to bind on oric in order to create a replication bubble. following the standard thermodynamic model of protein-dna binding, the probability of this event is dependent on the number of dnaa-atp molecules that are bound to the non-specific sites along the chromosome. these are low-affinity sites compared to titration sites, but the affinity is high enough so that the protein spends most of its time bound to the genome. these sequence-independent interactions are typical of dna-binding proteins. as a consequence, the simplifying assumption is usually made that the key molecular players are bound to the dna either specifically or non-specifically. simply stated, this is just an implementation of the known fact that dna-binding proteins, besides binding tightly to their target sequences, are generally "sticky" for dna, in a sequence-independent manner. this implies that the timing of dna replication initiation in the cell is determined by the ratio dnaa-atp to non-specific binding sites on dna in the cell and effectively results in the computation of the amount of dna per cell). note that in this case the volume of the cell is not taken into consideration since it is not the change in concentration of dnaa or of dna that determines the initiation time. the same will apply to the binding of dnaa and rna polymerase to the dnaa promoter. the second assumption is that at the time of initiation this ratio will be the same, independently of the growth rate and thus the number of origins. we ask how the parameters of this model must vary in order for this assumption to hold in the range of doubling times between and minutes. in the absence of autoregulation, the only factor that contributes to a decrease in the ratio of dnaa-atp to non-specific binding sites is the increase in dna after dna replication has begun. the complete model also includes the autoregulation of dnaa expression by dnaa binding to its own promoter and dnaa-atp transformation into dnaa-adp through the rida process. it is assumed that this ode description is applicable to a single average-cell on time scales shorter than the length of the cell cycle. this hypothesis could be challenged for the shortest observable doubling times, but the formulation of the model is dictated by maximizing simplicity. the values of the parameters are all taken or estimated from the available experimental measurements. they are shown in table together with the sources. the parameters are fixed with the values in the table for τ = mins and then are able to vary for the other values of τ in order to fix the ratio of dnaa-atp:chromosome length at the moment of initiation. †: chiara sag-gioro, anne olliver, bianca sclavi: multiple levels of regulation in the growth rate dependence of dnaa expression, submitted. formulation of the model timing of replication we take into consideration the situation where the cell cycle repeats itself identically i.e. balanced, exponential growth. following cooper and helmstetter, at a time c + d after the initiation time, the cell divides, i.e. that time must be an integer multiple of the doubling time. thus, if τ is the doubling time of the cell and x is the initiation time, then we must have x+c+d=τ, where n is the integer number of times that τ divides c + d. n can be viewed as the number of overlapping rounds of replication, and 2n is the number of origins. thus, this equation reflects the phenomenon of overlapping replication rounds. figure shows how x varies with the doubling time of the cell. defining y as the time at which the chromosome completes replication, we have y=τ-d. promoter term the activity of the dnaa promoter is the source term for dnaa-atp. we describe it by the standard thermodynamic model first used by shea and ackers. we denote the promoter term as q, the number of non-specific binding sites on the chromosome as nns and the number of rna polymerase molecules as p. furthermore, we denote the number of dnaa-atp molecules as a-. we then use the assumption that the number of non-specific binding sites is proportional to the length of dna in the cell, i.e. Λ = κnns. thus the expression for the rate of transcription at the dnaa promoter can be written as q=kaΘ1+c1Λp+c2a-p where ka is the basal rate of transcription of the dnaa promoter, Θ is the number of dnaa promoters in a given cell at a given time, and the remaining factor is the probability of rnap binding to the promoter. the parameters c and c depend on the binding energies Δεpd and Δεad of rnap and a- respectively to their promoter binding sites. the binding energies are determined from the ratio of specific vs non-specific binding affinities. c1=eΔεpd∕kbtκc2=e∕kbt. where the exponential terms are boltzmann weights. c = if the promoter is not autorepressed. a version of the promoter where dnaa binding to its sites is cooperative is described in additional file rida term this term reflects the number of dnaa-atp molecules that are converted to dnaa-adp molecules per unit time by the rida process. as discussed in the introduction, rida is a process that takes place at the replication forks during dna synthesis. we assume that the rate of conversion kr takes the same value at each replication fork. the number of pairs of replication forks at a given time, f, depends on which of x and y is larger. for x < y: f=2n-12⋅2n-12if0<t<xifx<t<yify<t<τ and for x > y: f=2n-12n-22if0<t<yify<t<xifx<t<τ and so the conversion from dnaa-atp to dnaa-adp takes place at a rate krf. this leads to the following differential equations for dnaa-atp and dnaa-adp ∂a-∂t=kaΘ1+c1Λp+c2a-p-krf ∂a+∂t=krf. term for the growth of the chromosome the growth of the chromosome is controlled by the replication forks. defining the rate of dna synthesis of each pair of replication forks as kΛ, we can write ∂Λ∂t=kΛf. assuming that kΛ is constant, and normalizing so that Λ = is the length of one full chromosome, we have kΛ = 1/c. main equation ∂r∂t=1ΛΘka1+c1Λp+c2Λpr-f. this equation describes the dynamics of the variable r=a-Λ which we suggest is a suitable candidate for the initiation potential since a specific number of dnaa-atp molecules is needed to be available to bind to the origin in order to induce dna melting, as described above. note that usually such dynamic equations are written in terms of volume, thinking of averages over cell populations on time-scales longer than a cell cycle. we assume that the background of genomic binding sites is the relevant variable in a shea-ackers type model, extending to a single cell cycle the approach normally used for longer time scales. we also assume that the volume can be treated as a weak perturbation, which we neglect here. in other words, the various molecules of interest are partitioned between the specific and non-specific binding sites on the chromosome. furthermore, the most important factor that determines the probability of binding to a given promoter is the absolute number of protein molecules relative to the absolute number of these binding sites, rather than the amount of protein per cell volume, which in comparison does not change significantly, and it can thus be neglected to a good approximation. thus, this assumption means that one need not track the volume of the cell, only the number of non-specific binding sites in the cell at a given time. this idea is discussed further in additional file section here we assume that the initiation potential, r, always reaches the same value at t = x independently of the growth rate and we ask how the parameters of this model must vary in order for this assumption to hold in the range of doubling times between and minutes. in the following, we will first establish that such a constant threshold cannot be obtained by a model with fixed parameters and then study the possible scenarios where different subsets of parameters are allowed to vary. main assumptions the model relies on the following further assumptions: all newly synthesised dnaa is immediately bound to atp, due to the relative abundance of atp in the cell compared to adp and the high affinity of dnaa for atp. dnaa-adp is only created by conversion from dnaa-atp by the rida process, when it is present. the probability of dnaa being bound to its sites at the origin or on the promoter, in the case of the presence of autorepression, is given by its thermodynamic equilibrium value. the same assumption holds for the binding of rna polymerase to a particular promoter. this means that we assume that the rate for transcription initiation is much slower than the rates for rnap binding and unbinding from the promoter. the rate of dnaa gene expression is proportional to the equilibrium probability that rnap is bound to the dnaa promoter. we do not consider translation directly and thus there is no time delay from transcription to protein production since the addition of this feature did not affect the result of the model. the number of non-specific binding sites on the dna in each cell for both rnap and dnaa is proportional to the total length of dna in the cell. the number of rnap molecules in the cell, p, grows exponentially from cell birth to cell division, corresponding to the hypothesis of constant concentration and exponentially growing cell size. a linear growth can also be used, however the dynamics of the model do not differ significantly between these two cases. numerical integration the non-linearity of the main equation necessitates the use of a numerical method of integration. we used a custom c++ implementation of the fourth-order runge-kutta method. the equation was integrated for values of the cell doubling time, τ, in the range mins ≤ τ ≤ mins. in order to test for the constant threshold condition, a transformation was performed by integrating the equation for τ = mins and using the value of r at initiation t = x as the imposed threshold for the other values of the doubling time. thus it was important to estimate, to a good degree of accuracy, values for the parameters at a doubling time of mins. these values appear in table the parameter values are either taken to be constant or are allowed to change and obtained as a consequence of the transformation. in all the scenarios, ka varies with growth rate. in scenario 1a, p changes as a function of growth rate with the values obtained from ref.. the dependence of the other parameters on growth rate is scenario specific.  a fixed set of parameters gives a varying initiation threshold with increasing growth rate we first describe the behaviour of the model with a fixed parameter set. the ratio r=a-Λ increases from the time of birth of the new cell. this can be interpreted as the accumulation of the'initiation potential'. at initiation, r peaks and then falls again due to the increase in the number of non-specific binding sites. when including the rida process, the total rida rate also increases following initiation, due to the higher number of active replication forks, contributing to the decrease in the initiation potential. however, the value of r at t = x varies for the different values of the doubling time t. thus, it appears that this model, with fixed parameters, cannot give a constant threshold that is reached at initiation in the range of growth rates considered here. this fact naturally leads to us consider a model in which the parameters are able to vary with growth rate. biologically, this is a natural requirement, as one may well expect from previous observations of the change in cellular components as a function of growth rate. a constant threshold condition implies alternative scenarios of growth-rate dependency in the circuit architecture the condition of a constant dnaa-atp/dna threshold at the time of initiation can be imposed by performing a mathematical transformation on the model and verifying the implications of this for the values of the parameters. the mathematical details of this transformation can be found in additional file section the transformation yields a fixed threshold in r that is reached at initiation. this can be seen upon comparing the plots in figure 3a and 3b. in the latter, the value of r is now the same at initiation for every cell doubling time. we estimate all initial values of the parameters from the literature and then we allow some of the parameters to change during the transformation. thus, the transformation procedure imposes a decision on which parameters to fix and which to allow to change. this determines different scenarios, as illustrated in figures and and also summarised in table. these plots explicitly show the required parameter changes at varying growth rates. in the first scenario, the rida rate, kr, is chosen to remain independent of growth rate, and the other parameters are allowed to change. the transformation fixes the scaling of ka with growth rate, and provides an equation that c and p must satisfy. since the constraints reduce to one equation with two unknowns, there is a set of possible solutions for c and p, and the scenario is underconstrained. one possibility is to fix one of the parameters to a particular trend, in turn determining the second parameter, determining different subscenarios. in the first of these sub-scenarios, 1a, the variation in the number of non-specifically bound rnap in the cell, p, with growth rate is fixed a priori by imposing the trend of a previous study, which partitions rnap into different classes, the rnap bound on promoters and non-specifically bound, and determines the dependence of this partition with growth rate,. imposing this trend determines the values for the binding affinity of rnap to the promoter, which must vary with growth rate. in sub-scenario 1b, the binding affinity of rnap to the dnaa promoter is chosen to remain independent of growth rate. fixing this parameter constrains the remaining equation, thereby imposing a particular dependence on growth rate of the number of non-specifically bound rnap, p, which turns out to be compatible with ref.. in scenario the binding affinity of both rnap to the dnaa promoter and of dnaa-atp to its repression sites are chosen to remain independent of growth rate, while the basal rate of transcription from the dnaa promoter, ka, the levels of free rnap, p, and the rida rate, kr are all permitted to vary from their original values. this scenario is also feasible in absence of autorepression. in brief, two possible scenarios can result in a constant initiation threshold for the model. in the first, the binding affinity of dnaa-atp to its repressor sites decreases with increasing cell doubling time and in the second the rida rate increases with cell doubling time. in both scenarios ka must decrease with increasing cell doubling time. it is important to note that ka, the basal dnaa transcription rate, needs to vary in all scenarios with growth rate. this term is independent from rnap availability and binding. it describes how quickly rnap moves through a gene when transcribing. the variation of this characteristic time with growth rate can be associated with variations in dna supercoling. in absence of rida or autorepression, the transformation can still be performed. however, it implies that the ratio c1/p should remain constant with growth rate. since p varies, this means that c would have to compensate exactly for the changes in p in order to keep a constant threshold. we have considered these further scenarios to be less probable, because they may result in a less robust control due to an unlikely fine-tuning of two parameters. the resulting scenarios are compatible with available knowledge on rnap availability and total dnaa expression given these scenarios, we have asked whether the predicted parameter variation with growth rate and the observables quantities produced by the model were compatible with the measurements and observations available in the literature. starting from the dependency of available rnap with growth rate, this is predicted and matched with available experimental data in ref.. scenario 1a assumes this dependency, and therefore automatically accounts for this observation. figure 4a demonstrates how also scenarios 1b and are broadly compatible with the results of klumpp and hwa since the average levels of non-specifically bound rna polymerase decrease with increasing doubling time in all cases. we now turn to the changes in measured expression of dnaa with growth rate. this can be measured by a reporter gene technique. figure 4b shows how the model predicts that the average expression of the dnaa gene should change with growth rate in all scenarios. this appears to be independent of the scenario chosen and is compatible with previous findings. we have also performed our own measurements, using a gfp reporter of the dnaa promoter encoded on a plasmid, see also additional file section 7), confirming this agreement. finally, at fixed growth rate, in order to determine whether this model reflects the main features of the regulatory network in the cell, we reproduced some of the experimental perturbations described in the literature. one of these experiments changed the rate of rida by changing the level of expression of the gene encoding for the hda protein, while others controlled the expression of dnaa independently of the dnaa promoter. we found that a constant threshold can be obtained upon a fold change in the rida rate. in order to determine how rida rate and autorepression strength influence the activity of dnaa-atp within a cell cycle in the model, we have monitored the amplitude of the oscillations of the ratio dnaa-atp:total genome length. the results show that increased autorepression contributes to a smaller amplitude of dnaa-atp within the cell cycle. a decrease in rida has the opposite effect, and is compensated for by an increase in autorepression. in either case the inclusion of these additional control factors appears to result in smaller amplitude in the oscillation of dnaa-atp and a smaller variation in the average amount of dnaa-atp per cell as the growth rate is varied. this may be advantageous for the use of dnaa as a transcription factor whose activity is responsive to changes in the replication status of the cell, via the rida and titration effects. finally, once we obtained a set of parameters that satisfies the constant threshold constraint, we modified the rida rate while leaving the other parameters unchanged, attempting to reproduce the effect of under or over expressing the hda protein, as in the experiments by riber et al.. in our model, this results in a change of the threshold value as a function of growth rate, and more significantly at slow growth. at faster rida rate, the threshold value is higher at slow growth, while the opposite is observed when rida rate is decreased. these compatibility tests give positive results, but do not allow us to distinguish between the two scenarios. we have explored the literature for tests of dependency of the rida rate with growth rate, and have found no evidence of this, which lead us to consider scenario where the rida rate is constant as the main one. model variants in order to gain confidence that the conclusion is not a consequence of the restricted set of biological ingredients included in this model, we considered some additional model variants, including some of the known factors that may influence the timing of replication initiation. delay in the synthesis of dnaa-atp we introduced a delay, representing the time necessary to obtain an active dnaa molecule from the binding of the rna polymerase to the dnaa promoter to the end of translation. this delay, however, does not produce a significant effect in imposing a dnaa-atp threshold at initiation, suggesting that translation delay might not have a predominant role in controlling the timing of initiation. cooperativity of autorepression cooperativity of autorepression affects the growth rate dependence of gene expression. additional file figure a6c shows by simulation that the presence of cooperativity alone cannot explain a constant threshold. thus it is still necessary to impose a constant threshold on the model using a transformation such as that described in additional file section however, in general one cannot make a translation and scaling and keep the promoter in the same form, and the transformation itself poses complex constraints on the possible regulations. thus, we decided to concentrate on the simpler non-cooperative model. biologically, of the two high-affinity dnaa sites found at the dnaa promoter, one matches exactly the consensus sequence, and has a higher affinity than the other, suggesting that at lower dnaa concentrations only one monomer could be bound. we have verified that the transformation is possible with a promoter that would fit this profile, i.e. of the form q=kaΘ1+c1Λp1+rk1+r2ωk1k21+rk where the parameter ω represents the cooperativity, and where k and k are the binding affinities of the two dnaa binding sites, multiplied by the proportionality constant between Λ and nns. this description introduces two new parameters - the cooperativity and the binding affinity of the second binding site. since these parameters appear in the equations describing the transformation that keeps the constant threshold, the scenarios become more underconstrained in a non-essential way. thus, while a more realistic promoter model could be useful for future descriptions, we decided not to pursue it here. the data locus and specific dnaa binding sites we considered the effect that the presence of the data locus has on the model by introducing a site on the chromosome that binds up to dnaa-atp molecules immediately after initiation has taken place. we included one data site for each origin in the cell. the results indicate that the data locus might indeed prevent further initiations in a given cell cycle, since it titrates large numbers of dnaa-atp molecules, effectively preventing them from binding at the origin. however, this variant of the model fails to achieve a constant dnaa-atp threshold at initiation at different growth rates if all parameters are kept constant. one can speculate that the large increase in data sites at faster growth rates would require a proportional increase in the rate of dnaa synthesis that cannot be solely provided by the increase in gene copy number. it must be also noted that the data locus appears to be unnecessary to prevent reinitiation and limit initiation of to once per cell cycle. additionally, we considered the effect of the reported ≈ binding sites distributed around the chromosome. these high-affinity sites would sequester dnaa in a similar way as the data locus, but proportionally to the genome amount, and thus effectively decrease its concentration. this is again insufficient to guarantee a constant initiation threshold and is equivalent to rescaling the rida rate. the eclipse no constraint for the eclipse period was included in this model. this proved not to be necessary, probably because there is no delay between the attainment of the threshold and the initiation of dna replication, thus avoiding the possibility that there is an'overshoot' of initiation potential before both rida and the increase of non-specific dna can begin. on the other hand, the gene copy number immediately doubles upon initiation of dna replication which results in a sudden increase in the number of dnaa-atp synthesized per unit time. however, this does not result in an increase in the initiation potential due to the corresponding increase in the number of replication forks and thus on the rida rate and the number of non-specific titrating sites. dnaa-atp recycling regions genomic recycling regions catalyzing the reconversion of dnaa-adp into dnaa-atp have been recently discovered. the quantitative contribution of this process is not clear, but within our framework it makes sense to ask how this would affect the initiation threshold, as this process is, roughly, a correction to rida. precisely, assuming dnaa-adp is not limiting, this would change the model by the substitution of kr with kr-ρΛ∕f where ρ is a fixed recycling rate. we verified that in the model this recycling term cannot by itself impose a constant threshold, while it can contribute to correcting quantitatively the effective rida rate. since none of these model variants qualitatively changes the behaviour of the model with respect to attaining a constant initiation threshold, they were not included in the minimal model formulation, in order to avoid confusion and proliferation of parameters. however, as shown by the variants explored above, the qualitative behaviour that the parameters of the model must vary with growth rate does not hold strictly for the minimal model only, but might be more general. discussion standard models of bacterial regulatory circuits were adapted to situations where the growth rate is fixed. the notion that these quantitative descriptions must account for bacterial physiology through the growth-rate dependent basic partitioning of the cell physico-chemical components is now entering the field of systems biology through a combination of new work and reconsideration of the classics. the dependency of the basic parameters on growth rate can produce notable effects on a genetic circuit, and complicates the standard descriptions. in our case, the task is more difficult, as the circuit under examination is active in determining some features of the bacterial physiology and not only affected by them. furthermore, on the technical level, one must produce a time-dependent description the expression of dnaa over cell-cycles of a range of durations. perhaps also for this reason, despite the fact that the regulation of dna replication has been a subject of intense study for over years, many questions remain open. given these obstacles, we have shown that, under a series of simplifying hypotheses, a consistent mean-field description for the dnaa/replication initation circuit is possible with varying growth rate. our description includes the processes that are believed to be most important for initiation of replication. in these respects, it is broadly compatible with previous modelling approaches. its originality lies in the minimality and in the attention given to growth-rate dependency. we focused on the minimal ingredients necessary in order for the basic tenet that the ratio dnaa-atp/dna attains a constant threshold at initiation to hold. the validity of this tenet is confirmed by the recent observations that initiation time is not affected by adding an extra origin on the chromosome and on the compensatory mutations emerging in hda mutants. we have defined the dna replication initiation potential, determining the timing of dna replication, as the dnaa-atp to dna ratio, r. molecular titration has been shown to result in ultrasensitive "all or none" responses, which further justifies using r as the threshold and could explain the synchrony of initiation in cells containing oric minichromosomes. we assume that its value at the time of initiation, r, is independent of the specific growth rate. the amount of dnaa-atp at the time of initiation thus needs to increase as a function of growth rate in order for r to remain constant as a function of doubling time, and we found that consequently, some of the model's parameter values must be allowed to vary. this assumption has not been verified directly. on the other hand, we feel that our point of view would still be useful in case of a growth-rate dependent r, as it is unlikely that this dependence would automatically match the dependence of all the other parameters. we have defined two main scenarios in which different subsets of the parameters are allowed to change. in scenario the rida rate, kr, is held constant as a function of growth rate, but the binding affinities of rnap and dnaa-atp to the dna need to vary with growth rate or allowing it to be free ). in scenario the binding constants are independent of growth rate but the rida rate, kr, must vary. we have verified that both scenarios are consistent with the eperimentally tested predictions of rnap availability with growth rate and with previous measurements and our own experimental evaluation of total dnaa expression, and also with a number of "in silico mutations" inspired by the available literature. thus, the scenarios appear as possibilities that are testable, but for the moment remain open. note that the property that the initiation threshold holds constant with respect to growth rate changes is not related to the specific set of parameters we used, or any set of parameters. our analysis shows that in general, for any fixed parameter set at a given growth rate, a transformation is necessary in order to keep the threshold constant while moving to another growth rate. in order to provide specific examples, we have produced plots in the style of those in figure with different curves corresponding to choosing different values of the initial input parameters. these demonstrate that the qualitative behaviour of the transformation is independent of these parameters. this exercise is also important to show that the parameter changes with growth rate are not numerically negligible for empirically plausible parameters, so that the question of keeping the initiation threshold constant is not purely academic. it is then interesting to ask which of these scenarios is more reasonable considering the known biological processes. we speculate that scenario is less likely, since until now there is no evidence pointing to a possible change in the intrinsic rida rate as a function of growth rate. the dnaa-related protein hda mediates this process. experiments with mutants over-and underexpressing hda, with corresponding increases and decreases in the rida rate, suggest a possible mechanism by which the kr term in the equations could vary by a growth rate-dependent expression of the hda gene. there may also be other, as yet unknown, factors that affect the growth rate dependence of the rida rate. alternatively, we can speculate that the decrease in the rate of rida with growth rate could be caused effectively by the action of the reverse process of dnaa-atp recycling by the recently discovered recycling regions. figure 5b shows that the rida rate should increase with cell cycle time and thus decrease with growth rate. this growth rate increase causes overlapping replication rounds, and thus higher chromosome copy number. since more recycling regions are present there is more recycling, i.e. a decrease in the effective rida rate, compatibility with the requirement imposed by our results. however, considering explicitly this model variant, we find that the balancing recycling cannot by itself impose a constant threshold. conversely in scenario the rida rate per replication fork is constant, and one has to rationalize the variation of the binding affinities. it seems possible that the binding affinities could change with growth rate through changes in supercoiling, in similar ways to those seen in figure and additional file figure a. the levels of average negative supercoiling are known to increase as the growth rate increases. however, it makes sense to challenge the validity of the basic assumption that the ratio of dnaa-atp to dna at the time of initiation is constant. this model assumes that the affinity for dnaa-atp binding to its own promoter can change with growth rate but its affinity for the origin does not. the first assumption mainly allows the model to change the magnitude of negative autoregulation as a function of growth rate, and it may indeed be explained by the changes in global cellular parameters such as negative supercoiling. we have considered how the activation threshold in the model would be affected if the binding affinity for dnaa to the origin would vary in the same way as its value at the dnaa promoter, required by scenario for a set of realistic parameters. we found that these changes in r are less than 10% over a wide range of growth rates, suggesting that this scenario might be robust. indeed, the observed threshold is certainly approximately constant when compared to the untransformed case i.e. the different values of the ratio at t = x shown in figure 3a. more generally, the initiation of dna replication has been significantly simplified in this model; all it requires is a specific amount of available dnaa-atp molecules. however we know that other factors, such as the binding of nucleoid proteins fis, ihf, h-ns and hu, may contribute to the formation of an open complex at the origin. on the other hand, other recent results have shown that at slower growth the cell contains a greater average amount of dnaa-atp per origin that results in initiation events that are independent of the novel synthesis of dnaa-atp. these results suggest that the regulation of the initiation process at the origin might indeed be dependent on the growth rate and that these changes still remain to be characterized quantitatively before they can be included in a theoretical model. interestingly, the basal rate of transcription of the dnaa gene, ka must vary in both scenarios. figure shows that ka decreases as the cells grow more slowly. this is what is expected from a promoter like the one of the dnaa gene that closely resembles ribosomal rna promoters. this family of promoters have a gc-rich sequence at the transcription initiation site called a discriminator region. this region renders the activity of the promoter sensitive to the degree of negative supercoiling, which activates transcription by enhancing dna melting, and leads to its inhibition by the accumulation of ppgpp at slower growth rates.  all things considered, we can say that perhaps our main result is that the determination of the timing of initiation by dnaa, besides relying on the known "architecture" comprising autorepression, rida and a number of other "dedicated" processes, can be understood only in its complex interplay with bacterial physiology nevertheless, it makes sense to ask whether this model allows us to elucidate some features of the reciprocal role of rida and dnaa autorepression, its two main ingredients. biologically, rida renders the control of dnaa-atp dependent upon ongoing dna replication, and thus results in an increase in dnaa-atp when replication forks are blocked. autorepression however probably plays a larger role in the absence of rida at slow growth, or in bacteria that do not have rida at all. the hda protein and thus the rida process seems to be quite specific to the fast-growing e. coli bacterium and its close relatives in the enterobacteriaceae, suggesting that in other bacterial species this level of regulation may not be required and is replaced instead by protein degradation, e.g. in caulobacter, or a high intrinsic atpase activity of the protein, as in mycobacterium tuberculosis. we have verified that the model can work in the absence of autorepression or rida, but the tuning of the parameters to achieve a constant threshold is more "difficult", in the sense that it requires more fine-tuning of the parameters, since the ratio c1/p should remain constant with growth rate. meaning that it is possible that a smaller range of growth rates would be accessible in these conditions. moreover, in the model, increasing autorepression or rida rate results in a smaller amplitude of the oscillations of the ratio dnaa-atp/dna during the cell cycle, and in a smaller variation in the average amount of dnaa-atp per cell as the growth rate is varied. this may be advantageous for the use of dnaa as a transcription factor which has to sense perturbations in the replication status of the cell at all growth rates. competing interests the authors declare that they have no competing interests. authors' contributions mcl, bs, and bb designed research. mg, bb, mcl, uf, and cs performed research. mg, bs and mcl wrote the paper. all authors read and approved the final manuscript. supplementary material additional file 1 additional text and figures. single pdf file containing the additional text and figures mentioned in the main text. click here for file acknowledgements we are grateful to matteo osella, rosalind allen, pietro cicuta, antonio celani, andrea sportiello, kunihiko kaneko and massimo vergassola for useful discussions and feedback. the authors acknowledge support from the human frontier science program organization and from epsrc.
in order to investigate the contributions of the different regulatory processes to the timing of initiation of dna replication at varying growth rates, we formulate a minimal quantitative model of the initiator circuit that includes the key ingredients known to regulate the activity of the dnaa protein. this model describes the average-cell oscillations in dnaa-atp/dna during the cell cycle, for varying growth rates. we isolate two main possible scenarios for this, depending on the roles of dnaa autoregulation and dnaa atp-hydrolysis regulatory process. we find that a quantitative description of replication initiation by dnaa must rely on the dependency of the basic parameters on growth rate, in order to account for the timing of initiation of dna replication at different cell doubling times. in escherichia coli, overlapping rounds of dna replication allow the bacteria to double in faster times than the time required to copy the genome. the precise timing of initiation of dna replication is determined by a regulatory circuit that depends on the binding of a critical number of atp-bound dnaa proteins at the origin of replication, resulting in the melting of the dna and the assembly of the replication complex. furthermore, both inactivation and autorepression reduce the amplitude of the cell-cycle oscillations of dnaa-atp/dna. we evaluate the conditions under which this ratio attains the same threshold value at the time of initiation, independently of the growth rate. in either case, the basal rate of gene expression needs to increase with the growth rate, in accordance with the known characteristics of the dnaa promoter. one possibility is that the basal rate of regulatory inactivation by atp hydrolysis must vary with growth rate. alternatively, some parameters defining promoter activity need to be a function of the growth rate.
 in escherichia coli, overlapping rounds of dna replication allow the bacteria to double in faster times than the time required to copy the genome. the precise timing of initiation of dna replication is determined by a regulatory circuit that depends on the binding of a critical number of atp-bound dnaa proteins at the origin of replication, resulting in the melting of the dna and the assembly of the replication complex. the synthesis of dnaa in the cell is controlled by a growth-rate dependent, negatively autoregulated gene found near the origin of replication. both the regulatory and initiation activity of dnaa depend on its nucleotide bound state and its availability.  in order to investigate the contributions of the different regulatory processes to the timing of initiation of dna replication at varying growth rates, we formulate a minimal quantitative model of the initiator circuit that includes the key ingredients known to regulate the activity of the dnaa protein. this model describes the average-cell oscillations in dnaa-atp/dna during the cell cycle, for varying growth rates. we evaluate the conditions under which this ratio attains the same threshold value at the time of initiation, independently of the growth rate.  we find that a quantitative description of replication initiation by dnaa must rely on the dependency of the basic parameters on growth rate, in order to account for the timing of initiation of dna replication at different cell doubling times. we isolate two main possible scenarios for this, depending on the roles of dnaa autoregulation and dnaa atp-hydrolysis regulatory process. one possibility is that the basal rate of regulatory inactivation by atp hydrolysis must vary with growth rate. alternatively, some parameters defining promoter activity need to be a function of the growth rate. in either case, the basal rate of gene expression needs to increase with the growth rate, in accordance with the known characteristics of the dnaa promoter. furthermore, both inactivation and autorepression reduce the amplitude of the cell-cycle oscillations of dnaa-atp/dna. 