A simple model for Lutz and Bujard's controllable promoters and its application for analyzing a simple genetic oscillator.

We develop an exact and flexible mathematical model for Lutz and Bujard's controllable promoters. It can be used as a building block for modeling genetic systems based on them. Special attention is paid to deduce all the model parameters from reported (in vitro) experimental data. We validate our model by comparing the regulatory ranges measured in vivo by Lutz and Bujard against the ranges predicted by the model, and which are calculated as the reporter activity obtained under inducing conditions divided by the activity measured under maximal repression. In particular, we verify Bond et al. assertion that the cooperativity between two lac operators can be assumed to be negligible when their central base pairs are separated by 22 or 32 bp [Gene repression by minimal lac loops in vivo, Nucleic Acids Res, 38 (2010) 8072-8082]. Moreover, we also find that the probability that two repressors LacI bind to these operators at the same time can be assumed to be negligible as well. We finally use the model for the promoter P(LlacO-1) to analyze a synthetic genetic oscillator recently build by Stricker et al. [A fast, robust and tunable synthetic gene oscillator, Nature, 456 (2008) 516-519].

Introduction

A defining goal of synthetic biology is to design and build genetic systems with specific functions that resemble those easily find nowadays in the electronic industry, such as genetic switches and oscillators. One of the main problems in the area is to design and characterize reliable synthetic promoters and operators that can control the expression of a particular gene via specific repressors, activators, or inducers. Lutz and Bujard [1] published in 1997 a seminal paper announcing the development of the synthetic controllable promoters P, P, P, P, and P in E.coli. For example, the promoters P and P have two sequences of the respective lac and tet operators (binding sites for the repressors LacI and TetR) inserted into regions adjacent to the -10 and -33 hexamers in the phage lambda promoter P; see the topography and sequences presented in Fig. 1. operators. So that, P and P are respectively repressed by LacI and TetR, and so induced by IPTG andanhydrotetracycline.

The EXPRESSYS company [2] is nowadays commercializing pZ vector systems that contain a module with one of the five Lutz and Bujard’s synthetic promoters followed by a ribosomal binding site and by either a multiple cloning site, the lacZ gene, or the luciferase gene. The pZ vector systems also have a second module with four possible origins of replication and a third module with one of five antibiotic resistance markers that carry their genuine promoters and ribosomal binding sites. The EXPRESSYS pZ vectors have been used by several authors [3-7] for in vivo experiments. In particular, Stricker et al. [4] used the pZ vectors for implementing in vivo a pair of genetic oscillators.

Considering the wide applications of Lutz and Bujard’s controllable promoters, it is important to explain in detail how they function; i.e., how polymerases and repressors are competing for occupying their respective binding sites in the promoter. In particular, it is important to explain why the in vivo regulatory ranges measured by Lutz and Bujard have the values reported in Table 1 of [1], because there are some apparent contradictions with several in vivo measurements done by Oehler et al. [8, 9] on the lac promoter. Recall that regulatory ranges are calculated as the reporter activity measured under inducing conditions divided by the activity measured under maximal repression. For example, the promoters P and P both contain a pair of lac operators and their in vivo regulatory ranges are 620 and 350, respectively. These values are at least one order of magnitude larger than the in vivo regulatory range (= 18) reported by Oehler, Eismann et al. in 1990 for a lac promoter with only one functional lac operator next to the transcriptional starting site; see configuration λEwt100 in Fig. 2 of [8]. Moreover, the regulatory ranges 620 and 350 are at least one order of magnitude smaller than the regulatory range (≥ 12 × 103) reported by Oehler, Amouyal et al. in 1994 for a synthetic lac promoter that includes two functional lac operators inside its structure; see Fig. 5 of [9].

Fig.2

Graphic representation of the different states O 0, , O and O 2 assumed by the promoter P in presence of the polymerase and the (dimeric) tetracycline repressor TetR. The transcriptional starting site is marked with the letter A, the -10 and -33 hexamers are indicated by small grey rectangles, the polymerase is represented by an orange kidney shaped figure, the tet operators are indicated by narrow red rectangles, and the repressors TetR are represented by big red ovals.

In this work we develop an exact and flexible mathematical model for Lutz and Bujard’s controllable promoters P, P, P, and P. Special attention is paid to deduce all the model parameters from reported (in vitro) experimental data. We propose a simple model (written in terms of chemical reactions) that explains how the polymerases and repressors compete for occupying their respective binding sites in the promoter. We also use this model for verifying Bond et al. assertion that the cooperativity between two lac operators can be assumed to be negligible when their central base pairs are separated by 22 or 32 bp [10]. Recall that two lac operators cooperate to produce a larger repression when a tetrameric repressor LacI binds simultaneously to them creating a repression loop [11], so that the cooperativity is proportional to the probability that a single repressor binds to two operators. Bond et al. and Müller et al. [10, 12] have found that, for lac operators whose central base pairs are separated by a distance between 49 and 98bp, the maxima and minima in the cooperativity are periodically distributed; i.e., the cooperativity is maximal when the separation between central base pairs is 59, 70, 81, or 92bp; while it is minimal (negligible) when the separation is around 52, 63, 74, 85, or 96bp. Moreover, Bond et al. also suggested that the non-coplanarity of the repressor LacI [11] repression loops impossible for distances shorter than 50bp, unless there are dramatic changes in the protein–DNA interface.

Finally, once we have a flexible mathematical model for Lutz and Bujard’s controllable promoters, we use it for analyzing the simplest synthetic genetic oscillator implemented in vivo by Stricker et al. in 2008 [4]. We analyze this oscillator because of its simplicity, it has only one feedback loop, so that the existence and frequency of oscillations are quite sensible to variations in the parameters. We find that the period of the simulated oscillations coincides with the period of the oscillations (≈35min) measured in vivo by Stricker et al. [4], so that we can take this coincidence as a second confirmation that our model for Lutz and Bujard’s controllable promoters is correct.

Methods and models

Consider the topography and sequences of the central parts of the lac promoter P and four of the Lutz and Bujard promoter/operators presented in Fig. 1.operators. The lac and tet operators are located besides (downstream) of the transcriptional starting site, besides (upstream) of the -33 hexamer, or between the -10 and -33 hexamers. Since the -10 and -33 hexamers are a principal landing mark for the polymerases in order to get attached to the promoter and to initiate the transcription, we use Chamberlin proposal that polymerases and repressors physically exclude each other from the promoter, and so we model this competition defining four basic states for the promoter: O0 represents the state when the promoter is free and ready to be occupied by a polymerase or a repressor; represents the state when a polymerase P is bound to the promoter but the DNA is still closed; O represents the state when one repressor R is bound to the respective operator (or to any one of the operators if there are two of them); O2 represents the state when there are two repressors bound to two operators. Notice that the state O2 only exists for the promoters P, P, and P, because only them contain two operators. For example, the different states found in promoter P are graphically represented in Fig. 2. states. The transitions between the states O0, and O are represented as chemical reactions:

where M represents a nascent mRNA molecule, the constant k is the promoter transcription initiation rate, and the parameters k, k, k, and k are the respective kinetic association and dissociation rate constants for the complexes repressor/operator O and polymerase/promoter . The parameter n denotes the number of active operators located in the core (central) part of the promoter and K2 is the association constant for the complex (2 repressors)/(2 operators) O2, so that n = 1 and K2 = 0 for the promoters P and P; while n = 2 for the other three promoters P, P, and P presented in Fig. 1. operators. Obviously, if there is only n = 1 operator in the promoter, then the state O2 does not exist and we may fix K2 = 0. Moreover, if there are n = 2 operators in the promoter, the resulting kinetic association constant nk must be twice as large as the original k, because the repressor has now two options (operators) to get bound to. The last reaction in (1) represents the transcription process, minus the initiation part; i.e.: promoter clearance, elongation, and termination. We must include it because there is a transition involved from the closed state into the free one O0. The reaction begins with a polymerase P recently bound to the operator (closed state ), the polymerase then opens the DNA strands, frees the promoter site (yielding the state O0), transcribes a nascent mRNA molecule M in the elongation process, and at the end detaches itself from the DNA.

We could not find experimental data on the value of K 2, and so we must use indirect information. Hillen et al. [13] have made several experiments weighting fragments of DNA that contain one or two tet operators, and in presence or absence of the repressor TetR. In Fig. 2 of their paper, Hillen et al. indicate that a fragment of DNA containing only one tet operator weights the equivalent of 95 base-pairs more when it is saturated with repressors TetR than when there are no repressors at all. The same experiment was repeated with a fragment of DNA that contains two tet operators. The central base pairs of these operators are separated by 29bp. Hillen et al. indicate that the fragment with two operators weights the equivalent of 174 base-pairs more when it is saturated with repressors TetR than when there are no repressors at all. Since 174 is approximately equal to 1.83 times 95, we can interpret the latter results as follows: When the fragment of DNA with two tet operators is saturated with repressors, the probability that it is actually bound by two repressors is 83% , while the probability that it is bound by only one repressor is 17% . Finally, since the central base pairs of the tet operators in the promoter P are separated by only 24bp, instead of 29bp, we make the following assumption: When the promoter P is mixed with repressors TetR, the probability that it is in the state O 2 is 75% , while the probability that it is in the state O is 25% .

We make a different assumption for the promoters P and P, which include a pair of lac operators inside their structure. The central base pairs of the operators are separated by less than 32.5bp according to Fig. 1. operators. Since the tetrameric repressor LacI has a molecular mass (= 154.52 kDa [11]) that is at least three times larger than the mass of the dimeric repressor TetR (= 46.71 kDa [14]), we shall assume that the state O 2 never happens for the promoters P and P; i.e., once a repressor LacI binds to one lac operator, it physically prevents a second repressor getting bound to the adjacent operator. The consistency of the previous assumptions is validated by comparing the regulatory ranges predicted by the resulting model against those measured in vivo by Lutz and Bujard. We rewrite below these assumptions to simplify the presentation.

Assumption 1. Let

This assumption means in particular that the probability P (O 2) is identically equal to zero for all the promoters analyzed in this paper, with the only exception of the promoter P. On the other hand, as we have said in the introduction, the tetrameric repressor LacI is able to bind simultaneously to two operators creating a repression loop. Hence, we should take into consideration the possibility that the lac operators inside the promoters P and P could cooperate to produce a larger repression. Bond et al. created in vivo a minimal system in E. coli for testing the shortest possible DNA repression loop between two lac promoters [10]. The authors suggested in particular that the non-coplanarity of the repressor LacI [11] makes repression loops impossible when the central base pairs of the promoters are separated by distances shorter than 50bp, unless there are dramatic changes in the protein–DNA interface. Since the central base pairs of the lac operators are separated by only 32bp (or 22bp) in the promoter P (or P), we made the followingassumption.

Assumption 2. The cooperativity is negligible between the

The assumptions made in the previous paragraphs automatically imply that we only need to consider the following states for the promoter: O 0, , O , and O 2. The last state O 2 is only needed for the promoter P. We must introduce some notation at this point. Let P (O 0), , and P (O ) be the probability that the promoter is free (state O 0), it forms the complex polymerase/promoter , and it forms the complex repressor/operator O , respectively. Moreover, we also need the probability P (O 2) that the promoter is at the state O 2. The state of equilibrium for the reactions (1) is easily given by the following identities,

Since , we can easily deduce the following equation after dividing by and applying the above identities and the considerations made in Assumption 1,

We can easily verify that the last term is really small (<10-3) with respect to the number 1, when we substitute the value of the parameters k , k , k , and [P] enlisted in (13), (14), (15), or (16). Notice that 1M is equal to 1.8 × 109 mpb (molecule count per average sized bacterium) according to the conversion calculated at the beginning of Appendix E. Whence, we can eliminate the last term in equation (3) without altering it, and so the probability that the promoter forms the complex polymerase/promoter is given by the formula

Finally, since the regulatory ranges are calculated as the reporter activity measured under inducing conditions divided by the activity measured under maximal repression, we can also calculate them as the probability obtained when [R] =0 divided by the probability obtained when the concentration [R] is maximal; i.e.,

It is known that the transcription process is a very complex one, which involves many consecutive chemical reactions; see for example the illustrative paper by deHaset et al. [15]. Nevertheless, even when the system of reactions (1) is a very simplistic way to represent the repression/transcription process, we assert that it is an effective and flexible model that can be used as a building block for modeling the genetic systems based on Lutz and Bujard’s controllable promoters. We compare in Section 3 the regulatory ranges measured in vivo for the promoters P, P, P, and P against those ranges obtained from equation (5) by substituting the value of the corresponding chemical kinetic constants measured in vitro by several authors. Since we obtain a good agreement between the measured and calculated ranges (the error is less than 7.7% in all cases), we can conclude that the system of reactions (1) and Assumptions 1 assumption and 2 are consistent with the experimental data.

Model for Stricker et al. genetic oscillator

As we have said in the previous paragraph, we obtain a good agreement between the measured and calculated regulatory ranges. However, since the regulatory ranges are calculated with data measured in steady state, we also want to verify that our model can make accurate predictions in a fully dynamical situation. To do so, we build a model for a synthetic genetic oscillator implemented in vivo by Stricker et al. [4]. The authors actually implemented two synthetic genetic oscillators using pZ vectors commercialized by the EXPRESSYS company. The topology of one oscillator has two feedback loops with different sign (that compete one against each other), while the topology of the second one is simpler and has only a single negative feedback loop [4]. We analyze the simpler oscillator, because its own simplicity means that the frequency and existence of oscillations are quite sensible to variations in the parameters. Elowitz and Leibler have already implemented in vivo a version of the repressilator (genetic oscillator with only one feedback loop) and observed that small variations in the parameters produce large variations in the dynamics [3]. Hence, an agreement between the period of the oscillations measured in vivo by Stricker et al. and the period calculated in silico with our model would be a strong confirmation that the system of reactions (1) can be used as a building block for modeling the genetic systems based on Lutz and Bujard controllablepromoters.

Stricker et al. implemented in vivo their simplest genetic oscillator by introducing two different EXPRESSYS pZ vectors into a strain of E.coli (denoted JS013). The first pZ vector contains a copy of the lacI structural genes controlled by the promoter P. Since the promoter is repressed by LacI and the lacI genes are transcribed and translated into the same repressor LacI, the production of LacI forms a simple negative feedback loop by repressing its own expression. The second pZ vector contains a yemGFP gene controlled by the promoter P, so that it works as a reporter for the activity of LacI. We begin our model by considering the system of reactions (1) that describe the dynamics (repression/transcription) of the P operator,

We do not include the state O 2, because it is only needed for the promoter P. The value of the parameters k , k , k , k , and k is deduced in Appendix C. We assume that the transcription is an instantaneous reaction, because the translation of the nascent mRNA molecule (M) can begin when the first two segments of the leader region have been transcribed and a ribosome can bind to them [11]. We add to our model the reactions corresponding to the translation of the lactose repressor monomer E; the dimerization of E into a functional molecule D that can bind to the lac operator; and the tetramerization into a complete (tetrameric) lactose repressor R that can bind to two differentoperators [11],

The parameter τ stands for the time delay involved in the translation process. The lactose repressor (LacI) is a tetrameric molecule, whose function as a repressor is blocked by the inducer 1-isopropyl-β-D-thiogalactoside (better known as IPTG). Properly speaking, the tetrameric lactose repressor R is a V-shaped molecule, whose branches are two dimers; each dimer D is able to bind independently to a lac operator; and molecules of IPTG can bind to a dimer D decreasing by one thousand its affinity for the lac operator itself [11, 16]. The different iterations between the IPTG (represented by I) and the molecules D and R can be represented by the following reactions,

where D is a dimer contaminated with IPTG and R (or R 2) is a tetramer with one branch (or both branches) also contaminated with IPTG. The contaminated dimers can participate in the tetramerization process, so that we must add the following two reactions,

Moreover, the tetrameric repressor R with only one branch contaminated with IPTG and the dimer D can both bind to one of the lac operators in P, rendering the transcription impossible, so that we add the following two reactions,

It is important to point out that the propensity k for the formation of the complexes O and is half of the propensity 2k for the formation of the original complex O given in (6). Finally, we complete our model for Stricker et al. oscillator by including the degradation of the mRNA (M) and the different proteins,

The value of the chemical kinetic parameters involved in the system of reactions (6) to (11) is deduced in Appendices C and E. We include below a list of these parameters. The standard units are minutes (min) and “molecule count per average sized bacterium” (mpb),

Finally, Stricker et al. have found that the in vivo implementation of their simplest genetic oscillator displays strong variations in the period of the oscillations [4]. The average period is approximately equal to 30 minutes and Stricker et al. also explained in the supplementary material of their article that the variations in the oscillations are a consequence of both the biochemical noise and the fact that their genetic oscillator has only one feedback loop. These results indicate that we must simulate the system of reactions (6) to (11) using an adequate stochastic algorithm, which allows us to include in our simulations the effects of both the biochemical noise and the delayed reaction in (7) that represents the translation. We choose the rejection method of Barrio et al., which is a modification of the Doob-Gillespie algorithm [17, 18] and which was proved by Cai to be an exact stochastic simulation algorithm (SSA) able to simulate the correct time evolution of a chemical reacting system that involves delayed reactions [19]. The pseudocode for the rejection method of Barrio et al. is presented in Appendix F. We perform one million simulations with this SSA; and to do so we consider thirteen chemical species: O 0, , O , O , , P, M, E, D, D , R, R , R 2, which interact through the chemical system of reactions (6) to (11).

Now then, we look for evidence of oscillations, and if it is the case, for the average period of such oscillations. Hence, we calculate the average power spectrum for the active repressor molecular counting, because any oscillation appears as a local maximum located at a frequency different from zero in the power spectrum. Moreover, the average period is given by the frequency (  ≠ 0) where the local maximum is located. We proceed as follows in order to analyze the system of reactions (6) to (11): We perform one million trajectories x  (t) with the rejection method of Barrio et al. described in Appendix F. Each trajectory x  (t) is performed from zero to 500 minutes and it measures the molecular counting of the active lactose repressor tetramer R at the time t. We have observed that the transitory state always happens between the first 100 minutes. Hence, we account only for the last 300 minutes from all the trajectories in order to eliminate completely the transitory state; and then we sample them recording a measurement every minute (we obviously get 300 samples from every trajectory). Finally, we calculate the fast Fourier transform of every sampled trajectory and work out the corresponding average power spectra. I.e., if c  (n) are the coefficients obtained from the fast Fourier transform of the sampled trajectory x  (t), the average power spectrum is given by the formula, where N = 106 is the number of performed trajectories. The results are presented in Section 4.

Results on regulatory ranges

We verify in this section that equation (5) correctly predicts the regulatory ranges measured in vivo for the promoters P, P, P, and P. The error is less than 7.7% in the four cases. Recall that the parameter n in (5) denotes the number of active operators located in the core part of the promoter, while the parameter β = 0 for all the promoters, with exception of P, for which we assume that β = 3.

Regulatory range of the lac promoter P

Oehler, Amouyal, et al. measured in vivo a regulatory range of 200 for the lac promoter with only one functional lac operator next to the transcriptional starting site (P) and a concentration of 50 tetrameric repressors (200 monomeric subunits) per cell [9, p. 3350]. We obtain the value 215.2 from equation (5) after substituting the value of the parameters below; see Appendix A.

As we said in the introduction, Oehler, Eismann et al. also measured in vivo a regulatory range of 18 for the promoter P in the E.coli strain CSH 9 rec A [8, p. 974]. We obtain the value 18.1 from equation (5) after substituting the value of the parameters above, with the exception of the number of active repressors per cell [R], for which we take 4 mpb instead of 50 mpb. The concentration [R] =4 mpb seems to be compatible with the assessments done by Gilbert et al. and Müller-Hill et al. [20, 21], who indicate that the repressor LacI is synthesized about 5 to 10 molecules each generation.

Regulatory range of promoter P

Lutz and Bujard measured in vivo a regulatory range of 350 for the promoter P, when controls the luciferase reporter gene in a pZ vector [1, p. 1205]. This promoter is build from the bacteriophage promoter P by adding two functional lac operators, one of them between the -10 and -33 hexamers. Lutz and Bujard induced a constitutive overproduction of the repressor LacI by placing a copy of the lacI gene under the control of the mutant constitutive laci promoter [21, p. 1260]. We obtain the value 372 from equation (5) after substituting the value of the parameters below; see Appendix B.

Lutz and Bujard measured in vivo a regulatory range of 620 for the promoter P, when controls the luciferase reporter gene in a Pz vector [1, p. 1205]. This promoter is build from the coliphage-lambda promoter P by adding two functional lac operators, one of them between the -10 and -33 hexamers. Lutz and Bujard induced a constitutive overproduction of the repressor LacI by placing a copy of the lacI gene under the control of the mutant constitutive laci promoter [21, p. 1260]. We obtain the value 595.9 from equation (5) after substituting the value of the parameters below; see Appendix C.

Regulatory range of the promoter P

Lutz and Bujard measured in vivo a regulatory range 2535 for the promoter P, when controls the luciferase reporter gene in a Pz vector [1, p. 1205]. This promoter is build from the coliphage-lambda promoter P by adding two functional tet operators, one of them between the -10 and -33 hexamers. Lutz and Bujard induced a constitutive overproduction of the repressor TetR by placing a copy of the tetR gene under the control of the constitutive promoter P . We obtain the value 2459.8 from equation (5) after substituting the value of the parameters below; see Appendix D.

Lutz and Bujard also measured in vivo the regulatory ranges 3670 (or 5050) for the same promoter P, when controls the luciferase reporter gene in a Pz vector that contains the origin of replication P15A (or pSC101*) instead of ColE1 [1]. We think that these repression ranges were obtained in bacteria that have 6 × 103 (or 4 × 103) active polymerases.

Results for Stricker et al. simple genetic oscillator

We verify in this section that the system of chemical reactions (6) to (11) is able to produce oscillations with an average period of approximately 30 minutes, when we use the parameters given in (12). A period of 30 minutes would be a good indication that reactions (6) to (11) correctly model the simplest genetic oscillator build by Stricker et al., because they have found in vivo oscillations with the same average period [4]. We look for oscillations by calculating the average power spectrum and observing if there are local maxima located at a frequency different from zero.

In Fig. 3 we present several average power spectra that we calculated for the molecular counting of the active lactose repressor tetramer R under different concentrations of IPTG. One can observe local maxima located at frequencies of approximately 0.5 and 0.6 miliHertz for IPTG’s concentrations of 100mpb and 150mpb, respectivelly. These local maxima correspond to oscillations with average periods of 33.3 and 27.8 minutes, and so we may conclude that there is a good agreement with the in vivo results obtained by Stricker et al. Recall that oscillations with a frequency of one milliHertz have a period of 16.67 minutes. It is also important to observe from Fig. 3 that the intensity of the oscillations diminish as the concentration of IPTG increases. This effect happens because the IPTG is indeed capturing the active lac repressor tetramers R and rendering them inactive.

On the other hand, recall that Stricker et al. implemented in vivo their oscillator by introducing two different pZ vector into a strain of E.coli. One pZ vector contains a copy of the lacI structural genes controlled by the promoter P; while the other contains a yemGFP gene controlled by the same kind of promoter. Thus, once we have verified that the system of reactions (6) to (11) correctly models the simplest genetic oscillator constructed by Stricker et al., we may wonder what would be the effect of substituting the promoter P by any other one.

Therefore, we now calculate the average power spectra for the molecular counting of the active lactose repressor tetramer R, when the promoter P has been substituted by P or P. To do so, we made again all the calculations indicated at the end of Subsection 1 model-Stricker, but we use the value of the parameters k , k , k , k , k and n indicated in (13) or (14) instead of those given in (12). In particular, one can verify that the average power spectra have the same shape and local maxima when the promoter P is substituted by P, so that the in vivo versions of the corresponding genetic oscillators should have similar behavior. Nevertheless, one also finds that the oscillations almost disappear behind the biochemical noise when the promoter P is substituted by P. In Fig. 4 we present four average power spectra calculated with the promoter P and different concentrations of IPTG. One can observe that the local maxima are not as evident as in Fig. 3-a, so that the oscillations are erratic and hidden by thenoise.

Conclusions

The positive results presented in previous Section 3 and 4 imply that the premises made in Assumptions 1 assumption and 2 are compatible with the experimental data, in particular, we can conclude that:

If the central base pairs of two lac operators are separated by less than 32.5bp, the probability that two repressors bind to them at the same time can be assumed to be negligible. I.e., once a repressor LacI binds to one of the lac operators, it physically prevents a second repressor getting bound to the adjacent one.

We verify Bond et al. assertion [10] that the cooperativity between two lac operators can be assume to be negligible when their central base pairs are separated by 22 or 32bp.

The repression and transcription initiation processes for Lutz and Bujard’s synthetic controllable promoters: P, P, P, and P can be modeled as the system of chemical reactions (1); and these reactions can be used as a building block for modeling more complex genetic systems based on Lutz and Bujard’s promoters.

The model of chemical reactions (1) correctly predicts, via equation (5), the in vivo regulatory ranges of the promoters: P, P, P, and P, when one substitutes the given parameters by the corresponding chemical kinetic constants measured in vitro by several authors. In particular, the probability that a polymerase forms a closed complex with the respective promoter can be expressed as the Hill function introduced in (4); i.e., where [R] and [P] are the number of active repressors and polymerases, respectively, and the value of m is obtained by substituting into equation (4) the parameters given in (13) to(16):

We deduce that the in vivo simplest genetic oscillator of Stricker et al. [4] must have a similar oscillatory behavior if the promoter P is substituted by P. Recall that Stricker et al. oscillator is based on a pZ vector, which contains a copy of the lacI structural genes controlled by the promoter P. Moreover, no oscillations must be observed if the promoter P is substituted by P. These results must be verified experimentally (in vivo).

As we said in the introduction, Oehler, Eismann et al. measured in vivo a regulatory range of 18 for the promoter P [8]. This range is one order of magnitude smaller than the ranges reported by Lutz and Bujard (350 and 620), because the concentration of repressor LacI is 4 mpb in Oehler et al. experiments, while the concentration is 90 mpb in Lutz and Bujardexperiments.

Oehler, Amouyal et al. also measured in vivo a regulatory range greater than 12 × 103 for a synthetic lac promoter that includes two functional lac operators inside its structure [9]. Oehler et al. explained that the magnitude of this range (≥12 × 103) is a consequence of the cooperativity existing between both lac operators, and so it is one order of magnitude larger than the ranges reported by Lutz and Bujard (350 and 620), because in the promoters P and P the lac operators do not cooperate to produce a larger repression.

Acknowledgments

We are quite grateful to the reviewers for their observations and suggestions, they have been quite useful for writing down the final version of our paper.

Appendix A: Parameters measured in vitro for the promoter P

Complex (repressor LacI)/operator formation rate constant k . Riggs et al. measured 7 (±0.9) ×109/M/s in vitro [22, p. 408]. We take k  = 7 ×109/M/s.

Complex (repressor LacI)/operator dissociation rate constant k . Riggs et al. measured 6.2 (±1.3) ×10-4/s in vitro [22, p. 404]. We take k  = 6.2 × 10-4/s.

Complex polymerase/(promoter P) formation rate constant k . Lanzer and Bujard state the rate 2 (±1) ×10/M/s [23, p. 8975]. We take k  = 2 ×106/M/s.

Complex polymerase/(promoter P) dissociation rate constant k . Meiklejohn and Gralla measured a half-life time of 3.5 minutes for the wild-type complex in presence of CRP-cAMP [24, p. 771]. We obtain .

Transcription initiation rate k . The polymerase opens the DNA in the second transcription step. Malan et al. measured an opening rate of 0.003/s for the lac wild-type promoter in presence of CRP-cAMP [25, p. 891], so we take k  ≈ 0.003/s.

Number of repressors LacI per cell [R]. Oehler et al. induced 200 monomeric subunits/cell [9, p. 3350], it corresponds to [R] =50 tetrameric repressors/cell.

Number of polymerases per cell [P]. Bremer and Dennis state there are between 1,500 and 11,400 polymerases/cell [26]. We take [P] =8300 polymerases/cell.

Appendix B: Parameters measured in vitro for the promoter P

Complex repressor/operator formation rate constant k  = 7 ×109/M/s; see Appendix A.

Complex repressor/operator dissociation rate constant k  = 6.2 × 10-4/s; see Appendix A.

Complex polymerase/(promoter P Lanzer and Bujard measured 6 (±0.5) ×107/M/s for the promoter P with a lac operator located between the -10 and -33 hexamers [23, p. 8975]. We take k  =  ×107/M/s.

Complex polymerase/(promoter P Brunner and Bujard measured k  = 9.2 × 10-4/s for the promoter P [27, p. 3141].

Transcription initiation rate k . The polymerase opens the DNA in the second transcription step. Dayton et al. measured a DNA opening rate of 0.09/s for the bacteriophage promoter P [28, p. 1618], so we take k ≈0.09/s.

Number of repressors LacI per cell [R]. Müller et al. state that the mutant promoter laci induces between 80 and 100 repressors/cell [21, p. 1260]. We take [R] =90 repressors/cell.

Number of polymerases per cell [P] =8300 polymerases/cell; see Appendix A.

Appendix C: Parameters measured in vitro for the promoter P

Complex repressor/operator formation rate constant k  = 7 ×109/M/s; see Appendix A.

Complex repressor/operator dissociation rate constant k  = 6.2 × 10-4/s; see Appendix A.

Complex polymerase/(promoter P Lanzer and Bujard measured 1.4 (±0.2) ×107/M/s for the promoter P with a lac operator located between the -10 and -33 hexamers [23, p. 8975]. We take k  = 1.4 × 107/M/s.

Complex polymerase/(promoter P Brunner and Bujard measured k  = 6.8 × 10-4/s for the promoter P [27, p. 3141].

Transcription initiation rate k . McClure observed a delay of 60s in the opening of complex polymerase/(promoter P ) [29, p. 5636]. We take .

Number of repressors LacI per cell [R] =90 repressors/cell; see Appendix B.

Number of polymerases per cell [P] =8300 polymerases/cell; see Appendix A.

Appendix D: Parameters measured in vitro for the promoter P

Complex (repressor (k /k ). Hillen et al. measured a value between 1012/M and 1013/M [13, p. 715]. We take k /k  = 5 ×1012/M.

Complex (repressor k . Hillen et al. measured a half-life of 12 min. [13, p. 715]. We obtain .

Complex polymerase/promoter formation rate constant k  = 1.4 × 107/M/s; see Appendix C.

Complex polymerase/promoter dissociation rate constant k  = 6.8 × 10-4/s; see Appendix C.

Transcription initiation rate ; see Appendix C.

Number of repressors TetR per cell [R]. Lutz and Bujard measured a ratio of 7/3 as the number of repressors TetR per cell divided by the number of repressors LacI [1]. Since there are 90 repressors LacI per cell according to Appendix B, we take [R] =210.

Number of polymerases per cell [P] =8300 polymerases/cell; see Appendix A.

Appendix E: Parameters involved in Stricker et al. simple oscillator

Some parameters were already given in Appendix C above; the remaining necessary parameters are estimated as follows:

E.coli volume. Volkmer and Heinemann measured the volume of a single E. coli under different growth conditions, they reported a variation in the volume from 1.6 to 4.4fl [30]. We take a mean volume of 3fl, so that we get the conversion

Number of lacI structural genes per cell There is only one lacI gene in the corresponding pZ vector. Stricker et al. (2008) state that there are between one and three pZ vectors per cell [4], so we take two structural lacI genes per cell.

Translation initiation rate for the lactose repressor monomer k . Chen et al. state the value k  = 10/min [31].

Dimerization and tetramerization formation rate constants k and k . We use the values k  = k  = 0.018/mpb/min given by Stricker et al. in [4] (suppl. material).

Complex (inducer IPTG)/(repressor LacI) formation rate constant k . Dunaway et al. measured in vitro the value k  = 4.6 × 104/M/s [32].

Dimerization and tetramerization dissociation rate constants k and k . They seem to be negligible in comparison to k  = k , so that we take k  = k  = 0/mpb/min.

Complex (inducer IPTG)/(repressor LacI) dissociation rate constant k . It seems to be negligible in comparison to k , so that we also take k  = 0/mpb/min.

Degradation rate constant for the lacI mRNA: γ . Varmus et al. measured a half-life time of approximately one minute [33, p.2264], so that .

Degradation rate constants for the repressor LacI: γ , γ , and γ . We assume that the monomeric E, dimeric D, and tetrameric R forms of the repressor LacI all have the same degradation rate. Chen et al. measured it, γ  = γ  = γ  = 0.1/min [31].

Time delay in the transcription of the lac mRNA molecule τ . We assume that τ  = 0, because the translation of the nascent lac mRNA molecule can begin when the first two segments of the leader region have been transcribed and a ribosome can bind to them [11].

Time delay in the translation of the lactose repressor monomer τ . The lactose repressor monomer is a protein of 360 amino acids residues [11, p. 523]. Bremer and Dennis states that the peptide chain elongation lies between 12 and 21 residues per second [26], so that .

Appendix F: Algorithms for simulating well-stirred chemical reacting systems

A well-stirred chemical reacting system is naturally described as a set of n ≥ 1 different chemical species {S 1, …, S } that interact through m ≥ 1 different chemical reactions. The dynamic state of such a system is represented by a (state) vector of natural numbers , so that the molecule count at time of the k-th chemical species (S ) is represented by the k-entry (X  (t)) of the vector X (t). The chemical reactions taking place in the reactor are described in the usual form:

where a (j, X) ≥0 denotes the propensity of the j-th chemical reaction, and is its corresponding stoichiometric vector. The Doob-Gillespie algorithm 17, 18] is an exact stochastic simulation algorithm (SSA) used for calculating the time evolution of the state vector X (t) under some assumptions: the chemical reactions in (17) are all instantaneous, no two chemical reactions can happen simultaneously, and the probability that the j-th reaction happens at the time interval [t, t + δ] is proportional to a (j, X (t)) δ. The pseudocode for the Doob-Gillespie algorithm is the following one:

Initialization. Set the time t = t 0 and the state vector X (t 0) = X 0 to their given initial values.

Calculate , where a  : = a (j, X (t)) is the propensity function for the j-th reaction in (17). Generate the waiting time from a standard uniform random variable u 1 in the interval [0, 1].

Generate a standard uniform random variable u 2 in the interval [0, 1] and choose the index μ = 1, 2, . . . , m such that .

Set the state vector X (t + τ) = X (t) + v [μ] and update the time t to t + τ. Go to the second step or else stop.

Since one can only analyze instantaneous reactions with the original Doob-Gillespie algorithm, it has been modified to produce several SSAs able to simulate the time evolution of a chemical reacting system that includes delayed reactions; i.e., chemical reactions like gene transcription and translation that take some time to deliver their products. Cai proved in [19] that the rejection method of Barrio et al. described below is an exact SSA able to include delayed reactions, and so we use it for calculating the time evolution of the state vector X (t) of the chemical reacting system (6) to (11). For simplicity, we adapt the algorithm for the case when there is only one non-consuming delayed reaction (indexed as the 0-th reaction)

By convention, in non-consuming reactions the reactants of an unfinished reaction can participate in new reactions, and so the corresponding molecular counts only change when the reaction is finished, like in gene transcription. The pseudo-code for the rejection method of Barrio et al. is the following one:

Initialization. Set the time t = t 0 and the state vector X (t 0) = X 0 to their given initial values.

Calculate , where a  : = a (j, X (t)) is the propensity function for the j-th reaction in (17) and (18). Generate the waiting time from a standard uniform random variable u 1 in the interval [0, 1]. If there is any delayed reaction (18) to finish in the time interval [t, t + τ), find the time t  ∈ [t, t + τ) of the first delayed reaction to finish, discard τ, set the state vector X (t ) = X (t) + v [0], update the time t = t , and repeat step 2. If there is no delayed reaction to finish in [t, t + τ), proceed to step 3.

Generate a standard uniform random variable u 2 in the interval [0, 1] and choose the index μ = 0, 1, . . . , m such that . If μ = 0, add the ending time t + τ 0 to the list of delayed reactions that are in process of finishing. If μ ≠ 0, set the state vector X (t + τ) = X (t) + v [μ].

Update the time t to t + τ. Go to the second step or else stop.

Notice that we need to include a first-in, first-out queue in this algorithm to save the ending times of the delayed chemical reactions that are in process of finishing.