Cycling expression and cooperative operator interaction in the trp operon of Escherichia coli

ARTICLE IN PRESS Journal of Theoretical Biology 263 (2010) 340–352

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Cycling expression and cooperative operator interaction in the trp operon of Escherichia coli Areli Herna´ndez-Valdez a, Moise´s Santilla´n b,c,, Eduardo S. Zeron d a

Universidad de Guadalajara, Centro Universitario de Ciencias de la Salud, Sierra Mojada 201, 44348 Guadalajara Jal, Me´xico ´n y de Estudios Avanzados del IPN, Unidad Monterrey, Vı´a del Conocimiento 201, Parque de Investigacio ´n e Innovacio ´n Tecnolo ´gica, Centro de Investigacio 66600 Apodaca NL, Me´xico c Centre for Nonlinear Dynamics in Physiology and Medicine, McGill University, 3655 Promenade Sir William Osler, Montreal, Canada QC H3G 1Y6 d ´ticas, Av. Instituto Polite´cnico Nacional 2508, 07360 Me´xico DF, Me ´n y de Estudios Avanzados del IPN, Departamento de Matema ´xico Centro de Investigacio b

a r t i c l e in f o

a b s t r a c t

Article history: Received 8 September 2009 Received in revised form 1 December 2009 Accepted 2 December 2009 Available online 22 December 2009

Oscillatory behaviour in the tryptophan operon of an Escherichia coli mutant strain lacking the enzymeinhibition regulatory mechanism has been observed by Bliss et al. but not confirmed by others. This behaviour could be important from the standpoint of synthetic biology, whose goals include the engineering of intracellular genetic oscillators. This work is devoted to investigating, from a mathematical modelling point of view, the possibility that the trp operon of the E. coli inhibition-free strain expresses cyclically. For that we extend a previously introduced model for the regulatory pathway of the tryptophan operon in Escherichia coli to account for the observed multiplicity and cooperativity of repressor binding sites. Thereafter we investigate the model dynamics using deterministic numeric solutions, stochastic simulations, and analytic studies. Our results suggest that a quasi-periodic behaviour could be observed in the trp operon expression level of single bacteria. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Systems biology Computational biology Nonlinear dynamics Stochastic simulation Time-delay differential equations

1. Introduction The tryptophan operon of Escherichia coli is one of the most extensively studied molecular regulatory systems. The detail with which the control mechanisms in this system are known makes it an excellent candidate to investigate its behaviour from a mathematical modelling standpoint. A number of models have been analysed in recent years to investigate different aspects of dynamics of this system (Santilla´n and Mackey, 2001; Venkatesh et al., 2004; Santilla´n and Zeron, 2004, 2006; Sima~ o et al., 2005; Bhartiya et al., 2006; Tabaka et al., 2008; Zeron, 2008; Nguyen and Kulasiri, 2008). In spite of that, there is still much to be learned from studying the trp operon. Bliss et al. (1982) carried out de-repression experiments with wild-type bacteria as well as with a mutant strain lacking the enzyme-inhibition regulatory mechanism. From their results the mutant bacterial strain showed clear sustained oscillations in mRNA, anthranilate synthase, and tryptophan intracellular

 Corresponding author at: Centro de Investigacio´n y de Estudios Avanzados del IPN, Unidad Monterrey, Vı´a del Conocimiento 201, Parque de Investigacio´n e Innovacio´n Tecnolo´gica, 66600 Apodaca NL, Me´xico. E-mail address: [email protected] (M. Santilla´n).

0022-5193/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2009.12.005

concentrations with a period close to 50 min.1 Later, Yanofsky and Horn (1994) also performed de-repression experiments with the same inhibition-free mutant bacterial strain. In contrast to Bliss et al., Yanofsky and Horn did not observe periodic oscillations. We also analysed with our previous model (Santilla´n and Zeron, 2004) the dynamic behaviour of the inhibition-free strain and were unable to reproduce the oscillatory behaviour reported by Bliss et al. Thus the possibility that the trp operon expression level in inhibition-free mutant E. coli bacteria oscillates remains unclear. On the other hand, it is known that three different operators (repressor binding sites) exist within the trpP promoter and that the first two of them interact cooperatively (Grillo et al., 1999; Jeeves et al., 1999; Brown et al., 1999). However, of all the models published so far, only that of Tabaka et al. (2008) takes this fact into account. In this last paper the authors focus on the dynamic behaviour of the operator–promoter DNA regulatory region and

1 They also developed a mathematical model which was able to reproduce their observed oscillations. Although time delays due to transcription and translation were taken into account, their model ignored transcription attenuation and, rather than modelling the repression mechanism in detail, considered a generic Hill-type equation with exponent equal to 4. Bliss et al. argued that the repressor is a homo-tetramer to justify this quite large exponent. However, it is now known that the repressor is a dimer and that no cooperativity exists between its subunits when they are bound by tryptophan. Hence, such a high Hill exponent cannot be justified from these considerations (Santilla´n, 2008).

ARTICLE IN PRESS ´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

conclude that it behaves as an efficient switch that responds accurately to fast consumption of tryptophan from the cell interior. In this paper we extend the trp operon model introduced by Santilla´n and Zeron (2004, 2006) to account for the observed operator multiplicity and cooperativity. The accuracy of the extended model is tested by simulating experiments that measure the dynamic response of the trp operon to nutritional shifts. Thereafter, we analyse the extended model response in a wild type strain and in a mutant strain that lacks the operon’s enzyme-inhibition regulatory mechanism. In particular, we are interested in investigating whether this mutation is capable of reproducing the oscillatory behaviour observed by Bliss et al. (1982). This is done using deterministic and stochastic numeric simulations, as well as analytically studying the model equations. If observed, the tryptophan operon cyclic behaviour could be important from the standpoint of synthetic biology whose goals include the engineering of genetic oscillators (Elowitz and Leibler, 2000; Wang et al., 2007; Li et al., 2007; Ninfa et al., 2007; Stricker et al., 2008). The paper is organized as follows. In Section 2 we introduce a detailed description of the tryptophan operon regulatory pathway, develop the mathematical model, and explain our parameter estimation procedure. In Section 3 the numerical methods used to solve the model time-delay differential equations and to carry out the stochastic simulations are explained. In Section 4 we present and discuss the results obtained from our numeric simulations and analytic studies. Finally, the conclusions are presented in Section 5.

2. Theory 2.1. The tryptophan operon The amino acid tryptophan can be synthesized by bacteria like E. coli through a series of catalysed reactions. The catalysing enzymes in E. coli are made up of the polypeptides encoded by the tryptophan operon genes: trpE, trpD, trpC, trpB, and trpA, which are transcribed from trpE to trpA. Transcription is initiated at promoter trpP, which is located upstream from gene trpE. The trp operon is regulated by three different negative-feedback mechanisms: repression, transcription attenuation, and enzyme inhibition. Below these regulatory mechanisms are briefly reviewed based on Yanofsky and Crawford (1987), Yanofsky (2000), Xie et al. (2003), Grillo et al. (1999), Jeeves et al. (1999), and Brown et al. (1999). The reader should consult Fig. 1 for a better understanding. Repression in the trp operon is mediated by three operators (O1, O2, and O3) overlapping with the operon promoter, trpP (see Fig. 1 A). When an active repressor is bound to an operator it blocks the binding of a mRNA polymerase to trpP and prevents transcription initiation. The trp repressor normally exists as a dimeric protein (called the trp aporepressor) and may or may not be complexed with tryptophan (Trp). Each portion of the trp aporepressor has a binding site for tryptophan. The trp aporepressor cannot bind the operators tightly when not complexed with tryptophan. However, if two tryptophan molecules bind to their respective binding sites the trp aporepressor is converted into a functional repressor. The resulting functional repressor complex can tightly bind the trp operators and so the synthesis of tryptophan catalysing enzymes is prevented. This sequence constitutes the repression negativefeedback mechanism: an increase in the concentration of tryptophan induces an increase in the concentration of the functional repressor, thus preventing the synthesis of tryptophan.

341

Transcription attenuation works by promoting an early termination of mRNA transcription, see Fig. 1B. The transcription starting site in the trp operon is separated from trpE by a leader region responsible for attenuation control. The transcript of this leader region consists of four segments (termed Segments 1, 2, 3, and 4) which can form three stable hairpin structures between consecutive segments. After the first two segments are transcribed they form a hairpin which stops transcription (cf. Fig. 1B-I). When a ribosome binds the nascent mRNA, it disrupts Hairpin 1:2 and transcription is re-initiated along with translation. Segment 1 contains two tryptophan codons in tandem. If there is scarcity of tryptophan, and thus of loaded tRNATrp , the ribosome stalls in the first segment. The development of Hairpin 2:3 (the antiterminator) is then facilitated, and transcription proceeds until the end (cf. Fig. 1B-II). However, if tryptophan is abundant, the ribosome rapidly finishes translation of Segments 1 and 2 and promotes the formation of a stable hairpin structure between Segments 3 and 4 (cf. Fig. 1B-III). mRNA polymerase molecules recognize this hairpin structure as a termination signal and transcription is prematurely terminated. Enzyme inhibition takes place through anthranilate synthase, the first enzyme to catalyse a reaction in the catalytic pathway that leads to the synthesis of tryptophan from chorismate. This enzyme is a hetero-tetramer consisting of two TrpE and two TrpD polypeptides. Anthranilate synthase is inhibited by tryptophan through negative-feedback. This feedback inhibition is achieved when the TrpE subunits in anthranilate synthase are individually bound by a tryptophan molecule, see Fig. 1C. Therefore, an excess of intracellular tryptophan inactivates most of the anthranilate synthase protein avoiding the production of more tryptophan. 2.2. Mathematical model A number of models for the tryptophan operon regulatory pathway have been published lately (Santilla´n and Mackey, 2001; Venkatesh et al., 2004; Santilla´n and Zeron, 2004, 2006; Sima~ o et al., 2005; Bhartiya et al., 2006; Tabaka et al., 2008; Zeron, 2008; Nguyen and Kulasiri, 2008). Most of them employ Hill-type functions to model the operon regulatory mechanisms. Only the models by Santilla´n et al. take into consideration all the involved chemical reactions to derive the regulatory functions. Hence, they are an optimal starting point to mechanistically account for the multiplicity and cooperativity reported by Grillo et al. (1999), Jeeves et al. (1999) and Brown et al. (1999). In the trp operon model developed by Santilla´n and Zeron (2004, 2006) the dynamics of the mRNA (M), anthranilate synthase (E), and tryptophan (T) concentrations are governed by the following differential delay equations: dM ¼ kM DRR ðTjttM ÞRA ðTjttM ÞðgM þ mÞM; dt

ð1Þ

dE 1 ¼ kE MjttE ðgE þ mÞE; dt 2

ð2Þ

dT T mT: ¼ kU ðtÞ þ kT ERI ðTÞr dt Kr þ T

ð3Þ

Here and thereafter, Xjtt denotes variable X delayed a time t, i.e. Xjtt ðtÞ ¼ XðttÞ. The delay tM is the time interval between transcription initiation and translation initiation, while tE is the time necessary to fully translate a TrpE polypeptide. Parameter D in Eq. (1) denotes the concentration of trp promoter copies. The factor 12 in Eq. (2) accounts for the fact that every anthranilate synthase protein contains two TrpE polypeptides. Parameter kM stands for the maximum rate of transcription initiation per promoter, kE is the maximum rate of translation initiation per mRNA, and kT is the

ARTICLE IN PRESS ´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

342

Fig. 1. Schematic representation of the three regulatory mechanisms found in the tryptophan operon: (A) repression, (B) transcription attenuation, and (C) enzyme inhibition. A glossary with the meaning of all the geometric forms in this figure is shown in panel D. See the main text for further explanation.

maximum rate of tryptophan production per enzyme E. Parameters gM and gE correspond to the M and E degradation rates. Parameter m is the bacterial growth rate. The term rT=ðKr þ TÞ represents the tryptophan consumption rate and kU ðtÞ stands for the influx of tryptophan from the environment. Finally, the functions RR ðTÞ, RA ðTÞ, and RI ðTÞ are all nonlinear monotone decreasing functions of T, respectively denoting the three different regulatory mechanisms present in the tryptophan operon: repression, transcription attenuation, and enzyme inhibition. Santilla´n and Zeron (2004, 2006), using chemical kinetics and making quasi-steady state assumptions for all fast processes, derived the following expressions for RR ðTÞ, RA ðTÞ, and RI ðTÞ:

RR ðTÞ ¼

P KP

 2 ; P RTot T 1þ þ 1 KP KR T þ KT

ð4Þ

T 1 þ 2a KG þ T RA ðTÞ ¼  2 ; T 1þ a KG þ T  RI ðTÞ ¼

KI T þKI

ð5Þ

2 :

ð6Þ

P represents the intracellular mRNA polymerase (mRNAP) concentration, RTot stands for the total repressor concentration, KP is the dissociation constant for the mRNAP–promoter complex formation reaction, KR is the dissociation rate for the repressor–operator complex formation reaction, KT is the dissociation constant for the reaction in which a tryptophan binds one of its corresponding binding sites in the aporepressor, a is a constant associated to the strength of transcription attenuation, KG is the dissociation rate for the tryptophan2tRNATrp complex formation reaction, and KI is the

ARTICLE IN PRESS ´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

3. Methods

Table 1 Parameter values as estimated in Appendix B. KR1 C0:625 mpb,

kcop C 11:13 KP C 2:0 mpb

P C 1500 mpb,

r C6:0  105 mpb min-1 m A ½0:01; 0:03 min-1

KI C 5000 mpb,

gE C0:01 min-1 , kM C5:1 min-1 ,

gM C 0:7 min-1 , Kr C 5000 mpb, R C 400 mpb,

We carried out both deterministic and stochastic simulations to analyse the model dynamic behaviour. The deterministic simulations were performed by numerically solving the set of time delay differential equations (1)–(3). To do this we used the fourth-order Runge–Kutta method implemented in xppaut (Ementrout, 2002). We also developed a stochastic simulation algorithm (SSA) to conduct the stochastic simulations. This SSA is a variation of Gillespie’s (1977) algorithm. Following Haseltine and Rawlings (2002) it partitions the chemical equations into slow and fast and makes an adiabatic approximation for all fast reactions. Moreover, it uses the ideas of Cai (2007) to account for the time delays. A more detailed description of the SSA implemented here can be found in Appendix C.3

KR2 C7:9 mpb

KR3 C100:0 mpb,

kE C30 min-1 KT C 20; 000 mpb D C 2 mpb a C 18:5 kT C 6:0  105 min

KG C 2500 mpb,

tM C 0:027 min-1 ,

tE C0:42 min

The meaning of mpb is ‘‘molecule count per average sized bacterium’’.

dissociation rate for the reaction in which a tryptophan molecule binds one of its binding sites in the anthranilate synthase enzyme. To take into account the existence of the three different repressor binding sites and the fact that the repressors that bind O1 and O2 do so cooperatively,2 we modified the function RR ðTÞ in Eq. (4) as follows:

RR ðTÞ ¼ 1þ

R2T KR1

! 1þ

R2T KR2

! 1þ

R2T KR3

P !KP þ

R22T KR1 KR2

1þ

! ; R2T P ðk 1Þ þ cop KP KR3

ð7Þ where R2T ¼ RTot



T T þKT

343

2

is the amount of active repressor, each KRi ði ¼ 1; 2; 3Þ is the dissociation constant for the R2T 2Oi complex formation reaction, and kcop 4 1 is a constant accounting for the cooperativity between Operators O1 and O2. The derivation of this regulatory function is explained in detail in Appendix A. Eqs. (1)–(3), together with the functions (5)–(7), constitute the complete set of delay differential equations modelling the regulatory pathway of the trp operon in E. coli which takes into consideration the existence of multiple operators within the promoter, as well as their cooperative behaviour.

2.3. Parameter estimation We paid special attention to the estimation of all the model parameters from available experimental data. The detailed procedure is described in Appendix B. The resulting parameter values are tabulated in Table 1.

2 Consider a receptor, A, and two ligands, B1 and B2 , that bind the receptor at two different sites. Let Ki and DGi , respectively, denote the dissociation constant and the Gibbs free energy of the reaction where ligand Bi binds the receptor independently of the other ligand. It is known that Ki ¼ expðDGi =RTÞ (Lehninger et al., 2008). If both ligands bind the receptor cooperatively, the free energy of the complex formed by the receptor and the two ligands is DG12 ¼ DG1 þ DG2 þ DGcop , where DGcop is due to the mutual interaction between ligands B1 and B2 . Finally, an effective dissociation constant can be defined for the overall reaction as K12 ¼ expðDG12 =RTÞ ¼ K1 K2 =kcop , with kcop ¼ expðDGcop =RTÞ. Since DGcop r 0 and the more negative this interaction energy the stronger the cooperativity between B1 and B2 , it follows that kcop Z 1 and that larger kcop values imply stronger cooperativity between both ligands.

4. Results 4.1. Deterministic simulations We started by analysing two different cases for the regulatory repression function RR ðTÞ. This was done to better understand the influence that the existence of multiple operators within the trp promoter and their cooperative interaction have on the operon dynamical behaviour. The first simplified situation corresponds to the existence of only one operator (Eq. (4)). In this case we took  Þ C 60, in agreement with the KP C 350:0 mpb so that RR ð0Þ=RR ðTmax experimental results of Yanofsky and Crawford (1987). In the second case we accounted for the existence of three operators overlapping the trp promoter, as well as of cooperativity between the first two of them. Thus we considered the repression function given by Eq. (7) with the parameter values in Table 1. Both functions are plotted in Fig. 2. We also fitted them to decreasing Hill-type functions of the form: HðTÞ ¼ RR ð0Þ

Kn ; Kn þ Tn

where RR ð0Þ is the value at T ¼ 0 of the corresponding repression regulatory function, and K and n are fitting parameters. The results of this fitting analysis reveal that multiple operators and cooperativity increases the half saturation constant K (from  2500 to  4500 mpb) and the exponent n (from  1:8 to  2:5). To test the model accuracy we simulated the de-repression experiments reported by Yanofsky and Horn (1994). In those experiments a bacterial culture was first grown in a tryptophan rich medium for a time long enough to ensure that the trp operon genes are silent. Afterwards, the bacteria were shifted to a tryptophan-free medium to induce expression of the operon genes. The anthranilate synthase enzymatic activity was then measured periodically, using colourimetry. We mimicked these de-repression experiments by numerically solving the model equations for 200 min, to guarantee that the system reaches a steady state, with a tryptophan influx rate of kU ¼ 1:0  106 mpb min1 , that is large enough to shut off the operon genes. Then we set t ¼ 0 and numerically solved the model equations with kU ðtÞ ¼ 0, for 0 rt r 150 min, using the final values of the previous solutions as the new initial conditions. These simulations were carried out using the two model variations described above. Since it is known that the bacterial growth rate varies greatly depending on the medium conditions we considered both the lower and upper limits reported for this parameter, mmin ¼ 0:01 min1 and mmax ¼ 0:03 min1 . The results are plotted 3 The stochastic simulation algorithm was implemented in python and the code is freely available at /http://web.me.com/moises.santillan/papers/trp3opSM. tgzS

ARTICLE IN PRESS 344

´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

and compared with the experimental results of Yanofsky and Horn in Fig. 3, under the assumption that the activity of anthranilate synthase is proportional to the active enzyme

Fig. 2. Plots of the repression regulatory function ðRR ðTÞÞ vs. T for the model variation that considers the three operators within the trp promotes and their cooperative interaction (solid line), as well as for the model variation accounting for only one operator (dashed line).

concentration, ERI ðTÞ. The simulated curves coincide with the experimental data points quite well in all cases, thus indicating no clear superiority between the two models within the context of this data. To analyse the influence of the time delays due to transcription and translation on the system dynamics we repeated the previous simulations by setting tM ¼ tE ¼ 0. The plots (not shown) corresponding to those in Fig. 3 reach the steady state at much the same time, but do so monotonically thus not showing the overshoot observed before the enzymatic activity reaches its steady state value. We conclude from this that the time delays play an important role in the operon de-repression dynamics, particularly during the first 30 min transient. As discussed in the Introduction, the possibility that the trp operon expression level in inhibition-free bacteria oscillates remains unclear and the question of whether the increased maximal slope of the repression regulatory function can explain an oscillatory gene expression in such mutant strain becomes interesting. To simulate the inhibition-free mutant we increased the value of the tryptophan-anthranilate synthase binding dissociation constant ðKI Þ and carried out several de-repression simulations with different values of this parameter. We also

Fig. 3. Plots of enzyme activity (normalized to one) vs. time obtained from simulations of the de-repression experiment carried out with the model variation that accounts for the existence of operator O1 only (Plots A and C) as well as with the model variation accounting for multiple operators and cooperativity (Plots B and D). The simulations -1 in Plots A and B were carried out considering the lower limit for the bacterial growth rate ðmmin C0:01 min Þ, while those in Plots C and D where calculated using the growth rate upper limit ðmmax C 0:03 min-1 Þ. In all plots the experimental points from a couple of experiments carried out by Yanofsky and Horn (1994) are also plotted with upward and downward pointing triangles.

ARTICLE IN PRESS ´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

345

Fig. 4. Bifurcation diagrams for the variables M, E, and T calculated with the minimum (A) and maximum (B) possible growth rate values. Here, the KI is the bifurcation parameter. Very large values of this parameter imply that tryptophan has low affinity for anthranilate synthase and so this enzyme cannot be inhibited. Solid lines denote stable steady states, dashed lines represent unstable steady states, and dash-dotted lines denote the upper and lower limits of stable limit cycles.

performed numerous simulations by eliminating the time delays and/or by considering the function RR ðTÞ given by Eq. (4) (which considers only one operator within the trp promoter) to investigate the possible role of time delays and operator multiplicity and cooperativity in the generation of oscillatory behaviour. Our results can be summarized as follows:

 The extended model predicts oscillatory behaviour for the







inhibition-free mutant strain during de-repression experiments. The increased maximal slope of the repression function (due to the three operators and their cooperative behaviour) as well as the time delays (due to transcription and translation) are both essential for these oscillations to appear. When either one or both of them are absent the model solutions show no oscillations. As seen in Fig. 4, KI can be considered as a bifurcation parameter. As its value is increased above a threshold level, the previously stable steady state becomes unstable and is replaced by an stable limit cycle. The bifurcation diagrams for M, E, and T when the growth rate takes its minimal value ðmmin C 0:01 min1 Þ are shown in Fig. 4A. Notice that the bifurcation occurs at KI  8:9  107 mpb, and the oscillation amplitude increases as KI increases. When KI -1 the mRNA and enzyme oscillation amplitudes approach 0.4 and 10 mpb, respectively. The oscillation period also increases as KI increases; it is about 120 min at the bifurcation point and tends to about 300 min in the limit KI -1. In Fig. 4B we show the bifurcation diagrams corresponding to the maximum value of the growth rate ðmmax C0:03 min1 Þ. Now the bifurcation takes place at KI  2:0  107 mpb and, as before, the oscillation amplitude increases with increases in KI . The amplitudes of the mRNA and enzyme oscillations, respectively, tend to 0.9 and 19 mpb in the limit KI -1. The oscillation period near the bifurcation point is less than 50 min and tends to about 140 min as KI tends to infinity. We numerically explored what other parameters influence the system oscillatory behaviour. We found that all the production

Fig. 5. Plots of mRNA (A), anthranilate synthase enzyme (B), and tryptophan (C) count per average bacterial cell vs. time calculated by simulating a de-repression experiment in a inhibition-free mutant strain. For this simulation we used m ¼ 0:03 min-1 and KI ¼ 5:0  109 mpb.

rates (kM , kE , kT ), as well as all the degradation or consumption rates (gM , gE , and r) affect the period and amplitude of oscillations and change the KI bifurcation value. Decreasing kM , kE , or kT and/or increasing gM , gE , or r have the same effect as increasing the bacterial growth rate. That is, the oscillation period decreases, the oscillation amplitude increases, and the KI bifurcation value decreases.

ARTICLE IN PRESS 346

´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

The oscillatory behaviour predicted by the present model is illustrated in Fig. 5. There we show M, E, and T vs. t corresponding to the simulation of a de-repression experiment carried out with m ¼ 0:03 min1 and KI ¼ 5:0  109 mpb. The period and amplitude of the oscillations depend on the actual parameter values, as described above, but the shape of the different curves agree with the experimental results of Bliss et al. (1982). One feature worthy of note is that the oscillation amplitude is quite small, and this is always observed whenever oscillations occur. The amplitude of mRNA oscillations is less than 1 mpb at its maximal value, which raises the question of whether these oscillations would be masked by random fluctuations once the so-called biochemical noise is taken into consideration. 4.2. Stochastic simulations The molecule count for the species involved in the tryptophan operon regulatory pathway is quite low in some cases: each cell typically has a few or none mRNA molecules and only a few dozen enzymes. This makes the mass action law, on which the deterministic model depends, unreliable. Rather than following the smooth trajectories predicted by the deterministic model the system variables would show large fluctuations also known as biochemical noise. The smaller the number of

molecules, the larger the corresponding fluctuations. See Raser and O’shea (2005) and Maheshri and O’shea (2007) for extensive discussions about biochemical noise in the context of gene expression. We carried out stochastic simulations of the model dynamic behaviour to take this biochemical noise into account. For this we implemented a variation of Gillespie’s algorithm, described in Section 3 and in Appendix C, which is capable of dealing with the time delays due to transcription and translation and simulates stochastically the evolution of the mRNA, enzyme, and tryptophan molecule counts. To understand the dynamic influence of the time delays due to transcription and translation we carried out stochastic simulations including and ignoring them. In Fig. 6 we show the T vs. t plots corresponding to simulations in a single wild-type bacterium as well as the averages over 100 simulations (in an equal number of independent wild-type bacteria). Notice that the single-bacterium simulations are quite similar. Both predict a mean tryptophan count around 20,000 molecules with a standard deviation of about 7000 molecules. The averaged plots reveal that, as expected, the fluctuation amplitudes are largely reduced in both cases. Indeed, the standard deviations decrease to about 800 molecules. On the

Fig. 6. Stochastic simulation of the tryptophan-count time evolution within a single wild-type bacterium (panels A and B) and average of the number of tryptophan molecules in 100 independent wild-type bacteria (panels C and D). These simulations were performed with the model variation that accounts for the existence of three operators within the promoter and of cooperativity among them. Furthermore, in the simulations corresponding to panels A and C we ignored the time delays due to transcription and translation, while we took these delays into consideration in the simulations of panels B and D.

ARTICLE IN PRESS ´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

other hand, a large initial overshoot is observed in the T vs. t plot when the time delays are included but not when they are ignored. This overshoot survives the averaging process because a quite-similar initial fluctuation is present in all single-cell simulations. In contrast, all the other fluctuations in single-cell simulations cancel because they have random amplitudes and occur at random times. Bhartiya et al. (2006) analysed the importance of this overshoot from an evolutionary standpoint and concluded that it allows bacteria to respond faster to nutritional shifts. We also performed stochastic simulations mimicking a mutant bacterial strain that lacks the enzyme inhibition regulatory mechanism. This was done by increasing the value of parameter KI in Eq. (6) up to 4:0  109 mpb. The results of these simulations are shown in Fig. 7. A comparison of the corresponding plots in Figs. 6 and 7 demonstrates that the tryptophan-count fluctuation amplitude increases enormously when enzyme inhibition is eliminated. Therefore, enzyme inhibition plays a major role in the reduction of biochemical noise. When the results of 100 independent simulations are averaged we observe, once more, a large initial overshoot in the T vs. t plot corresponding to the simulation with delays. The conclusion in this case is the same as in the simulations

347

with wild-type bacteria including delays (Fig. 6). Namely all single-cell simulations show the first pronounced fluctuation because of the time delays, and so it survives the averaging process. Observe that the plots in Figs. 7A and B are remarkably different. Although large fluctuations on the tryptophan molecule count can be noticed in both cases, the fluctuation amplitude is larger and apparently more regular in Fig. 7B. Note that the duration of fluctuations in Fig. 7A seems to be more disperse than those in Fig. 7B. To explore this fact further we calculated the duration of each fluctuation as the time between two consecutive local minima, computed the corresponding histograms, and plotted them in Fig. 8. We can see there that the two histograms differ greatly. On the one hand, the histogram in Fig. 8A (which corresponds to the simulation with no delays) is maximal at the origin and decreases monotonically until it vanishes at about 80 min. Conversely, the histogram in Fig. 8B (corresponding to the model with delays) is concave, maximizes at about 60 min, and vanishes at about 100 min. The fact that the histogram in Fig. 8B possesses a single maximum indicates the existence of a preferred fluctuation duration, though with high variability around it. We call this behaviour quasi-periodic.

Fig. 7. Stochastic simulation of the tryptophan-count time evolution within a single mutant bacterium lacking enzyme inhibition (panels A and B) and average of the number of tryptophan molecules in 100 independent mutant bacteria (panels C and D). These simulations were performed with the model variation that accounts for the existence of three operators within the promotes and of cooperativity among them. The value of parameter KI was increased up to 4:0  109 mpb to mimic the lack of enzyme-inhibition. Furthermore, in the simulations corresponding to panels A and C we ignored the time delays due to transcription and translation, while we took these delays into consideration in the simulations of panels B and D.

ARTICLE IN PRESS ´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

348

increasing functions of T. Moreover, after solving for T in (10), the solution TðEÞ is a growing function of E. The following equations are satisfied in the steady state: M ¼

E ¼

kM

gM þ m

DRR ðTðE ÞÞRA ðTðE ÞÞ;

kE M ; 2ðgE þ mÞ

ð11Þ

where M and E are, respectively, the M and E steady state values. Since both RR ðTÞ and RA ðTÞ are decreasing functions of T and TðEÞ is an increasing function of E, the equation system (11) has a unique solution ðM ; E Þ with both M  and E real and positive. The stability of the steady state is determined by the roots of the following characteristic equation (MacDonald, 1989): ðA þ lÞðB þ lÞ ¼ expðtlÞ; C

ð12Þ

with t ¼ tM þ tE , A ¼ gM þ m, B ¼ gE þ m, and   kM kE D dðRR RA Þ dT  C¼ : 2 dT TðE Þ dE E From their definition A; B 4 0 while C o 0. In case t ¼ 0 the characteristic equation (12) reduces to ðA þ lÞðB þ lÞ ¼ C:

ð13Þ

The left-hand side of this equation is a convex parabola with roots

l ¼ A and B. It follows from this and the fact that C o 0 that

Fig. 8. Histograms for the time periods between two consecutive minima of the tryptophan molecule count calculated from the two stochastic simulations shown in Fig. 7. Recall that those simulations correspond to bacterial strains lacking enzyme inhibition (KI ¼ 4:0  109 mpb). In the simulation corresponding to Panel A (that shown in Fig. 7 A) the time delays due to transcription and translation were ignored, while they were taken into account in the simulation corresponding to panel B (that shown in Fig. 7 B).

either (13) has two negative real roots or it has two conjugate complex roots with negative real part. In conclusion, the system steady state is always stable in the absence of time delays, regardless of parameter values (Strogatz, 1994). To test whether time delays can induce sustained oscillations, we assume that l ¼ io. Thus, when the left-hand side of (12) is plotted on the complex plane, the result is the lower branch of a parabola that opens to the right and has its vertex on the real axis at AB=C o0. On the other hand, the plot of the right-hand side of (12) is a circle with unitary radius. If both plots intersect, then the assumption that l ¼ io is correct and oscillatory solutions are possible. Otherwise, l is never purely imaginary and oscillatory solutions shall never be observed. A simple condition for the parabola and the circle to intersect is that AB ojCj. We evaluated constants A, B, and C and plotted the above referred parabola and unit circle for both the wild type bacterial

4.3. On the origin of oscillations Below we study the origin and behaviour of oscillations in the inhibition-free mutant strain of E. coli through a local stability analysis and further numeric simulations. Consider a system with no tryptophan uptake from the environment ðkU ¼ 0Þ. It has been proven that the dynamics of variable T are much faster than those of M and E (Santilla´n and Zeron, 2004, 2006). Therefore, the equation system (1)–(3) can be reduced via a quasi-steady state approximation to dM ¼ kM DRR ðTðEjttM ÞÞRA ðTðEjttM ÞÞðgM þ mÞM; dt

ð8Þ

dE 1 ¼ kE MjttE ðgE þ mÞE; dt 2

ð9Þ

where TðEÞ is the solution of kT ERI ðTÞ ¼ r

T þ mT: Kr þ T

ð10Þ

Eq. (10) has only one positive real solution because RI ðTÞ is a decreasing function of T while both terms in its right hand side are

Fig. 9. Unitary circle (solid line) and curves, respectively, determined by the rightand left-hand sides of Eq. (12) for the wild-type (dashed line) and the inhibitionfree mutant (dot-dashed line) bacterial strains.

ARTICLE IN PRESS ´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

strain (KI ¼ 5:0  103 mpb) and for a mutant strain lacking enzyme inhibition (KI ¼ 4:0  109 mpb) in Fig. 9. Notice that AB4 jCj, and so the parabola and the circle do not intersect for the wild type strain. This further means that time delays cannot destabilize the steady state to generate sustained oscillations. Conversely, AB o jCj for the inhibition-free mutant strain and so the circle and the parabola intersect, meaning that time delays can give rise to cycling behaviour. This change happens because jCj increases as KI increases, as can be tested by differentiating (10), solving for dT=dE, and substituting into the definition for C. The above analysis explains why time delays are necessary for sustained oscillations to emerge, and how this happens when KI is increased to simulate an inhibition-free bacterial strain. Let us now analyse the dynamics of oscillations and their behaviour in the presence to biochemical noise. When it cycles, the trp operon of E. coli behaves as a relaxation oscillator (Andronov et al., 1966). If anthranilate synthase concentration is high, the tryptophan level is high as well because enzymes constantly produce tryptophan due to lack of enzyme inhibition. Moreover, mRNA synthesis is effectively null since the dynamics of M are substantially faster than those of E. These facts imply on their own that intracellular enzyme concentration slowly decays with a rate depending on the bacterial growth rate and the enzymatic degradation rate. The system state does not stop at the steady state, because it is unstable, and E keeps decreasing. When the level of E becomes very low, the bacteria reach a state of tryptophan starvation. As a matter of fact, when T decreases below a given threshold the transcription initiation probability suddenly increases due to the increased sigmoidicity of function RR ðTÞ. Concomitantly, a few mRNA molecules are synthesized in a short period of time, and a few hundred enzymes are assembled when such mRNA molecules are translated. Thus, bacteria return to the state of high enzyme concentration and the process is repeated all over again. The oscillation amplitude predicted by the deterministic model is fairly small in all cases. Interestingly, when the system inherent biochemical noise is taken into account, the amplitude of oscillations increases by at least one order of magnitude. To analyse this behaviour recall that the amplitude of mRNA oscillations predicted by the deterministic model is less than one molecule (for instance it is 0.3 in the simulation of Fig. 5A). This occurs because the deterministic model relies upon the mass action law and the meaningful variables there are the chemicalspecies concentrations rather than the corresponding molecule counts. Contrarily, the stochastic model directly accounts for molecule numbers. Therefore, while the deterministic model predicts that a fraction of mRNA molecules are assembled once the transcription initiation probability gets close to unity, the stochastic model predicts the production of an entire number of them. This explains why oscillations, when they occur, have much lager amplitudes in the stochastic simulations than in the deterministic ones. On the other hand, when a given oscillation is ending the Trp concentration reaches very low values and this increases the transcription initiation probability. However, the exact time at which the next transcription is to be initiated (marking the beginning of the next oscillation) is a random variable due to the system stochastic nature. This last fact explains why the duration of fluctuations is only similar but not identical in the stochastic simulations.

5. Conclusions In this work we have extended and improved a mathematical model for the tryptophan operon previously introduced by

349

Santilla´n and Zeron (2004, 2006). This improvement consisted in modifying the repression regulatory function RR ðTÞ to account for the existence of multiple operators within the promoter region and of cooperativity between two of them (Grillo et al., 1999). We have also re-estimated some parameters using more accurate experimental data. The model accuracy was tested by simulating de-repression experiments carried out by Yanofsky and Horn (1994). The obtained results demonstrate that the improved model reproduces the experimental results as well as the original model, Fig. 3. To investigate the dynamic influence of time delays we repeated the above simulations after setting tM ¼ 0 min and tE ¼ 0 min. A comparison of the simulations with and without delays (results not shown) reveals that delays make the intracellular anthranilate synthase activity recover faster after a nutritional shift. In fact it experiences an overshoot before reaching the steady state value (Fig. 3). The existence of this overshoot was confirmed by stochastic simulations performed to study the influence of biochemical noise (Fig. 6). Bhartiya et al. (2006) assert that this initial overshoot gives E. coli an evolutionary advantage because it reaches the Trp necessary to grow and reproduce more rapidly. These facts prove that the time delays due to transcription and translation are very important dynamically speaking, regardless of their length. What is important in this respect is the place they occupy in the operon regulatory pathway. The overshoot discussed in the previous paragraph was not observed in the experiments by Yanofsky and Horn (1994) because of their time resolution. Our results suggest that this overshoot is not an artefact and that it may provide E. coli with an evolutionary advantage. Therefore, experimentally testing its existence becomes important. Below we examine an experimental setup that could be implemented to tackle this and other questions. As discussed in Section 4, there is a controversy regarding the oscillatory behaviour of an E. coli mutant strain that lacks enzyme inhibition. We mimicked this mutant strain by increasing the value of parameter KI and simulated de-repression experiments with it. Our objective was to test whether the improved model can reproduce the oscillatory behaviour originally reported by Bliss et al. (1982) but not confirmed by other studies (Yanofsky and Horn, 1994). According to the results of our deterministic simulations the trp operon is capable of showing sustained oscillations, but only when the time delays and the operator multiplicity and cooperativity are included. Nonetheless, the amplitude of oscillations is extremely low in all cases (for mRNA it is about 1 mpb at most). Interestingly, the removal of a negative feedback loop gives rise to an oscillatory behaviour. This apparently contradicts the common knowledge that negative feedback loops are necessary for oscillations to appear (Nova´k and Tyson, 2008). Nonetheless, there is no contradiction because, according to a linear stability analysis we performed on the model, repression and transcription attenuation are directly responsible (together with the time delays due to transcription and translation) for the predicted cycling behaviour. We carried out stochastic simulations to investigate the influence of biochemical noise on the oscillatory behaviour predicted by the deterministic simulations. Our results demonstrate that quasi-periodic oscillations could be observed in single bacteria but not in bacterial populations because of the oscillation-period variability induced by biochemical noise (Fig. 7). In this case, time delays are also essential for this quasi-periodic behaviour to emerge. Interestingly, biochemical noise not only does not make oscillations fade out, as one might expect, but

ARTICLE IN PRESS 350

´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

actually enhances them. The price to pay is the variability on the oscillation period mentioned above. Let us go back to the controversy regarding the oscillatory behaviour of E. coli’s inhibition-free strain. Yanofsky and Horn (1994) did their de-repression experiments by taking periodic samples from a bacterial culture and using colourimetry to determine anthranilate synthase activity. This technique involves numerous cells in every measurement. Therefore, it is not surprising according to our results that Yanofsky and Horn could not observe oscillations, even if they existed. On the other hand, Bliss et al. (1982) experiments also involved several cells in single measurements and they reported observing oscillations in the trpoperon expression level of an inhibition-free E. coli strain. We have not enough information to explain why Bliss et al. detected oscillations. All we could do is to speculate that perhaps something in their experimental set-up synchronized all bacteria in the culture. In any case, our model results suggest that quasi-periodic oscillations should be observed in individual inhibition-free bacteria. If experimentally confirmed, this oscillatory behaviour could be important in the context of synthetic biology, whose goals include the development of genetic oscillators (Elowitz and Leibler, 2000; Wang et al., 2007; Li et al., 2007; Ninfa et al., 2007; Stricker et al., 2008). In particular, Stricker et al. (2008) recently engineered a genetic relaxation oscillator in E. coli with tuneable oscillatory periods. The trp operon of an inhibition-free E. coli strain could be another example for this type of oscillators; its frequency could be tuned by simply changing the growing medium conditions and so the bacterial growth rate. Our results further agree with the assertion of Nova´k and Tyson (2008) that the necessary condiments for sustained oscillations are negative feedback regulation, increased nonlinearity of the regulatory functions, and a delayed response. As we have seen, the multiplicity of operators and their cooperative behaviour increase the nonlinearity of the repression regulatory function. Hence, operator multiplicity and cooperativity (or any other feature that increases the nonlinearity of the repression or attenuation regulatory functions) are necessary conditions, together with the time delays due to transcription and translation, for the predicted oscillatory behaviour. Interestingly, Bacillus subtilis has a tryptophan operon which shares the same regulatory pathway architecture with that of E. coli, although the specific mechanisms differ between both bacteria (Gollnick et al., 2005). Thus, It would be interesting to investigate whether these two systems present similar dynamics. Among them, oscillatory expression of the operon genes in inhibition-less bacterial strains. The initial overshot of the operon expression level after a nutritional shift, and the oscillatory behaviour of inhibition-less mutant bacteria, are predictions of the present work which demand experimental proof. For that, and experimental setup with high temporal resolution and the capability of measuring individual cells is necessary. This could be achieved by inserting a reporter gene (like gfp), controlled by the trp promoter, into the bacterial chromosome and continuously recording fluorescence in individual bacteria.

Acknowledgements This research was partially supported by Consejo Nacional de Ciencia y Tecnologı´a (CONACyT, MEXICO) under Grant: 55228. We are grateful to Prof. Michael C. Mackey for his critical reading of the manuscript and his useful comments and criticisms. We also thank the anonymous referee for his recommendations, which helped us to improve the paper.

Appendix A. Modelling of the repression function As seen in Fig. 1A one mRNA polymerase molecule (P) and up to three active repressors ðR2T Þ can bind the trp DNA regulatory region in the repression regulatory mechanism through the following reactions: KR1

R2T þ O!O1R ;

KR2

KR3

R2T þ O!O2R ;

KR2 =kcop

KR3

R2T þ O1R ! O12 2R ; KR3

123 R2T þ O12 2R !O3R ;

R2T þ O!O3R ;

R2T þ O1R !O13 2R ;

KR3

R2T þ O2R !O23 2R ;

KP

P þO!OP :

ð14Þ

In these reactions O, OiR , Oij2R , O123 3R , and OP , respectively denote the states where all operators are free, where only operator Oi is bound by an active repressor, where operators Oi and Oj are both bound by active repressors, where all three operators are repressor-bounded, and where a polymerase is bound to the promoter. Furthermore, KRi and KP are the dissociation constants for the Oi : R2T and promoter–polymerase complex formation reactions, respectively. Finally kcop 4 1 denotes the cooperativity between operators O1 and O2. The equilibrium equations for the reactions in Eq. (14) are R2T O ¼ KR1 O1R ;

R2T O ¼ KR2 O2R ;

R2T O1R ¼ KR2 O12 2R =kcop ;

R2T O ¼ KR3 O3R ;

R2T O1R ¼ KR3 O13 2R ;

R2T O2R ¼ KR3 O23 2R ;

3 123 R2T O12 2R ¼ KR O3R :

ð15Þ

Here and thereafter we employ the same symbol to denote a chemical species and its concentration. The conservation equation for the total trp DNA-regulatory-region concentration is 13 23 123 Oþ O1R þO2R þ O3R þO12 2R þ O2R þ O2R þ O3R þ OP ¼ OTot :

ð16Þ

Clearly, the reaction list in (14) does not contain all the possible reactions between active repressors, polymerases, and the trp DNA regulatory region. However, the equilibrium equations corresponding to the missing reactions can be recovered as linear combinations of the equilibrium equations in Eq. (15). In that sense, Eqs. (15) and (16) form a complete system of algebraic equations for all the possible operator-state concentrations. By solving for OP in Eqs. (15) and (16) we obtain, after some algebra, the following expression for the probability that a regulatory region is bound by a polymerase: OP RR ðTÞ ¼ ¼ OTot

1þ

R2T KR1

! 1þ

R2T KR2

! 1þ

R2T KR3

P !KP þ

R22T KR1 KR2

1þ

! : R2T P ðkcop 1Þ þ 3 KP KR

Appendix B. Parameter estimation From Grillo et al. (1999) an active repressor molecule R2T can bind three different operator sites (O1, O2, and O3) all of which overlap the trp promoter. Additionally, when two repressors are bound to O1 and O2 they interact in such a way that the absolute value of the binding energy is larger than the sum of the absolute values of the binding energies of single repressors separately bound to O1 and O2. Grillo et al. also report several measurements from which the following dissociation constants can be estimated: KR1 C 0:625 mpb; KR2 C 7:9 mpb;

ARTICLE IN PRESS ´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

KR3 C 100 mpb: The abbreviation mpb stands for ‘‘molecule count per averagesize bacterium’’. Each KRi denotes the dissociation constant for the reaction in which a repressor molecule binds the operator Oi alone. Grillo et al. also measured the binding energy when operators O1 and O2 are simultaneously bound by repressor molecules, and the cooperativity constant can be estimated from their experiments: kcop C11:13: Since there are about 1500 mRNA polymerase molecules per bacterium (Bremmer and Dennis, 1996), and the operon activity decreases 60 times due to repression (Yanofsky and Crawford,  Þ C 60 and the parameters P and KP are thus 1987), RR ð0Þ=RR ðTmax estimated as P C 1500 mpb; KP C2:0 mpb:  is the maximum steady-state intracellular concentration Tmax of tryptophan. According to Bliss et al. (1982)  C20; 000 mpb: Tmax

Caligiuri and Bauerle (1991) measured the apparent dissociation constant for anthranilate synthase inhibition by tryptophan. According to them half of the anthranilate synthase enzymes are inhibited when there are about 20,000 tryptophan molecules within a bacterial cell. Then from Eq. (6) KI C 5000 mpb: To estimate the maximum tryptophan consumption rate r note that: the average number of proteins per E. coli cell is about 2:6  106 , the average number of amino acids per protein is 360, and the abundance of tryptophan amino acids in proteins is around 1.53%. This means that there should be on the order of 1:43  106 trp amino acids assembled into the bacterium proteins.4 Since a bacterium can reproduce in 24 min under optimal conditions (Bremmer and Dennis, 1996) and the number of proteins has to be doubled before that happens, then the maximum rate at which tryptophan is consumed to synthesize protein is

r C6:0  105 mpb min1 : Bliss et al. (1982) and Sinha (1988) claim that the anthranilate synthase degradation rate is so low that it is very hard to measure in a bacterium lifetime. Hence, an upper bond for this parameter is

gE C 0:01 min1 : For the following parameters we used the estimates in (Santilla´n and Zeron, 2004) kM C 5:1 min1 ;

kE C30 min1 ;

gM C 0:7 min1 ;

KT C 20; 000 mpb;

Kr C5000 mpb;

D C2 mpb;

R C400 mpb;

a C 18:5;

KG C 2500 mpb: E. coli can reproduce every 24–60 min depending on the conditions of the growing medium (Bremmer and Dennis, 1996). Consequently, the bacterial growth rate can vary in the range

m A ½0:01; 0:03 min1 : 4

http://redpoll.pharmacy.ualberta.ca/CCDB/cgi-bin/STAT_NEW.cgi

351

With all the above parameters estimated, the value of kT (the rate of tryptophan production per anthranilate synthase enzyme) was calculated so that the average steady-state tryptophan concentration is about 20,000 molecules per bacterium when the culture is growing in a minimal media (Bliss et al., 1982): kT C6  105 min1 : To estimate the delay between transcription and translation initiation consider that just after a polymerase has transcribed the first two segments of the leader region it stops and waits for a ribosome to start translation. Since the added length of these two segments is 81 base pairs (Yanofsky, 2004) and the polymerase transcription rate is about 50 base pairs per second (Bremmer and Dennis, 1996) then

tM ¼

81 bp C0:027 min: 50 bp=s

Anthranilate synthase is a hetero-tetramer made up of two TrpD and two TrpE subunits. The genes trpD and trpE are, respectively, 514 and 531 codons long5 and they are translated separately. From these pieces of information and the fact that the average translation rate per ribosome is 21 codons per second (Bremmer and Dennis, 1996), the time delay due to translation is estimated as

tE ¼

531 codons C0:42 min: 21 codons=s

Appendix C. Stochastic simulation algorithm A stochastic simulation algorithm (SSA) was implemented to study the influence of biochemical noise on the operon dynamic behaviour. Given the large number of reactions involved in the trp operon regulatory pathway and the fact that most of them have comparatively large propensities, implementing the original Gillespie algorithm is not feasible. It would be too costly computationally speaking. Hence, a simplification is necessary in the reaction scheme. Following Haseltine and Rawlings (2002) the chemical reactions in the system were partitioned into slow (transcription, translation, and mRNA and enzyme degradation) and fast reactions (all the others). An adiabatic approximation was then assumed for all fast reactions. After doing that we ended with a reduced system involving the reactions tabulated in Table C1, together with their corresponding effective propensities. Of the reactions in Table C1, transcription and translation do not deliver their products immediately but take times tM and tE to do so, respectively. Recently, Cai (2007) introduces an algorithm to simulate the stochastic evolution of systems involving reactions with delays. Although faster stochastic simulation algorithms capable of dealing with delayed reactions have been proposed, the algorithm of Cai is exact in the same sense than that of Gillespie. This is why we decided to use it in our simulations. Cai’s SSA can be summarized as follows: 1. Initialize the system state, i.e. set the mRNA and TrpE polypeptide counts to their initial values. 2. Initialize the stack that will contain all unfinished reactions and the waiting times until they are finished. Set it empty. 3. Stochastically calculate the waiting time until the beginning of the next reaction, t. This computation is quite sophisticated because it involves various iterative steps during which all programmed reactions with a waiting time shorter than t are 5

http://www.ncbi.nlm.nih.gov/sites/entrez/

ARTICLE IN PRESS ´ndez-Valdez et al. / Journal of Theoretical Biology 263 (2010) 340–352 A. Herna

352

Table C1 Reduced reaction scheme for the trp operon regulatory pathway, after eliminating all fast variables via an adiabatic approximation. Reaction

Effective propensity

Transcription mRNA degradation Translation TrpE polypeptide degradation

kM DRR ðTÞRA ðTÞ ðgM þ mÞM kE M=2 ðgE þ mÞE

The value of T used in the calculation of RR ðTÞ and RA ðTÞ is obtained by numerically solving the following equation, given the value of E: kT ERI ðTÞ ¼ rT=ðKr þ TÞ þ mT.

4. 5. 6.

7.

8.

finished, the system state is concomitantly updated, and the reactions propensities are recalculated. See Cai (2007) for details. Update of the simulation time: t ¼ t þ t. Update the waiting times for all reactions programmed in the stack, i.e. subtract time t from every waiting time. Calculate which is the next reaction to begin. This reaction is calculated stochastically considering that the probability of each reaction is proportional to its propensity. If the reaction resulting from step 5 delivers its products immediately, update the system state accordingly. Otherwise, program it in the stack, together with its corresponding waiting time. Go to step 3 or finish the program.

References _ Fishwick, W., 1966. Theory of Oscillators. Andronov, A.A., Vitt, A.A., Khaı˘kin, S.E., Dover, New York. Bhartiya, S., Chaudhary, N., Venkatesh, K.V., Doyle, F.J., 2006. Multiple feedback loop design in the tryptophan regulatory network of Escherichia coli suggests a paradigm for robust regulation of processes in series. J. R. Soc. Interface 3 (8), 383–391. Bliss, R.D., Painter, P.R., Marr, A.G., 1982. Role of feedback inhibition in stabilizing the classical operon. J. Theor. Biol. 97, 177–193. Bremmer, H., Dennis, P.P., 1996. Modulation of chemical composition and other parameters of the cell by growth rate. In: Neidhart, F.C., Curtis, R., Ingraham, J.L., Lin, E.C.C., Low, K.B., Magasanik, B., Reznikoff, W.S., Riley, M., Schaechter, M., Umbarger, H.E. (Eds.), Escherichia coli and Salmonella thyphymurium: Cellular and Molecular Biology, vol. 2. American Society for Microbiology, Washington, DC, pp. 1553–1569. Brown, M.P., Grillo, A.O., Boyer, M., Royer, C.A., 1999. Probing the role of water in the tryptophan repressor–operator complex. Protein Sci. 8 (6), 1276–1285. Cai, X., 2007. Exact stochastic simulation of coupled chemical reactions with delays. J. Chem. Phys. 126, 124108. Caligiuri, M.G., Bauerle, R., 1991. Identification of amino acid residues involved in feedback regulation of the anthranilate synthase complex from Salmonella typhimurium. J. Biol. Chem. 266, 8328–8335. Elowitz, M.B., Leibler, S., 2000. A synthetic oscillatory network of transcriptional regulators. Nature 403 (6767), 335–338. Ementrout, B., 2002. Simulating, Analyzing, and Animating Dynamical Systems: A Guide to xppaut for Researchers and Students. Society for Industrial Mathematics, Philadelphia, USA. Gillespie, D.T., 1977. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81 (25), 2340–2361. Gollnick, P., Babitzke, P., Antson, A., Yanofsky, C., 2005. Complexity in regulation of tryptophan biosynthesis in bacillus subtilis. Annu. Rev. Genet. 39, 47–68.

Grillo, A.O., Brown, M.P., Royer, C.A., 1999. Probing the physical basis for trp repressor–operator recognition. J. Mol. Biol. 287, 539–554. Haseltine, E.L., Rawlings, J.B., 2002. Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics. J. Chem. Phys. 117, 6959–6969. Jeeves, M., Evans, P.D., Parslow, R.A., Jaseja, M., Hyde, E.I., 1999. Studies of the escherichia coli trp repressor binding to its five operators and to variant operator sequences. Eur. J. Biochem. 265 (3), 919–928. Lehninger, A.L., Nelson, D.L., Cox, M.M., 2008. Lehninger Principles of Biochemistry, fifth ed. W.H. Freemn, New York URL: /http://www.loc.gov/catdir/enhance ments/fy0835/2007941224-d.htmlS. Li, C., Chen, L., Aihara, K., 2007. Stochastic synchronization of genetic oscillator networks. BMC Syst. Biol. 1, 6. MacDonald, N., 1989. Biological Delay Systems: Linear Stability Theory. University Press, Cambridge. Maheshri, N., O’shea, E.K., 2007. Living with noisy genes: how cells function reliably with inherent variability in gene expression. Annu. Rev. Biophys. Biomol. Struct. 36, 413–434. Nguyen, L.K., Kulasiri, D., 2008. On multiple regulatory mechanisms in the tryptophan operon system in Escherichia coli: in silico study of perturbation dynamics. In Silico Biol. 8 (5–6), 485–510. Ninfa, A.J., Selinsky, S., Perry, N., Atkins, S., Xiu Song, Q., Mayo, A., Arps, D., Woolf, P., Atkinson, M.R., 2007. Using two-component systems and other bacterial regulatory factors for the fabrication of synthetic genetic devices. Methods Enzymol. 422, 488–512. Nova´k, B., Tyson, J.J., 2008. Design principles of biochemical oscillators. Nat. Rev. Mol. Cell. Biol. 9 (12), 981–991. Raser, J.M., O’shea, E.K., 2005. Noise in gene expression: origins, consequences, and control. Science 309, 2010–2013. Santilla´n, M., 2008. On the use of the Hill functions in mathematical models of gene regulatory networks. Math. Model. Nat. Phenom. 3 (2), 85–97. Santilla´n, M., Mackey, M.C., 2001. Dynamic regulation of the tryptophan operon: a modeling study and comparison with experimental data. Proc. Natl. Acad. Sci. 98, 1364–1369. Santilla´n, M., Zeron, E.S., 2004. Dynamic influence of feedback enzyme inhibition and transcription attenuation on the tryptophan operon response to nutritional shifts. J. Theor. Biol. 231, 287–298. Santilla´n, M., Zeron, E.S., 2006. Analytical study of the multiplicity of regulatory mechanisms in the tryptophan operon. Bull. Math. Biol. 68, 343–359. Sima~ o, E., Remy, E., Thieffry, D., Chaouiya, C., 2005. Qualitative modelling of regulated metabolic pathways: application to the tryptophan biosynthesis in E. coli. Bioinformatics 21 (suppl. 2), 190–196. Sinha, S., 1988. Theoretical study of tryptophan operon: applications in microbial technology. Biotechnol. Bioeng. 31, 117–124. Stricker, J., Cookson, S., Bennett, M.R., Mather, W.H., Tsimring, L.S., Hasty, J., 2008. A fast, robust and tunable synthetic gene oscillator. Nature 456 (7221), 516–519. Strogatz, S.H., 1994. Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering. Addison-Wesley, Reading, MA. Tabaka, M., Cybulski, O., Holyst, R., 2008. Accurate genetic switch in Escherichia coli: novel mechanism of regulation by co-repressor. J. Mol. Biol. 377 (4), 1002–1014. Venkatesh, K.V., Bhartiya, S., Ruhela, A., 2004. Multiple feedback loops are key to a robust dynamic performance of tryptophan regulation in Escherichia coli. FEBS Lett. 563 (1–3), 234–240. Wang, Z., Hou, Z., Xin, H., Zhang, Z., 2007. Engineered internal noise stochastic resonartor in gene network: a model study. Biophys. Chem. 125 (2–3), 281–285. Xie, G., Keyhani, N.O., Bonner, C.A., Jensen, R.A., 2003. Ancient origin of the tryptophan operon and the dynamics of evolutionary change. Microbiol. Mol. Biol. Rev. 67, 303–342. Yanofsky, C., 2000. Transcription attenuation, once viewed as a novel regulatory strategy. J. Bacteriol. 182, 1–8. Yanofsky, C., 2004. The different roles of tryptophan transfer RNA in regulating trp operon expression in E. coli versus B. subtilis. Trends Genet. 20, 367–374. Yanofsky, C., Crawford, I.P., 1987. The tryptophan operon. In: Neidhart, F.C., Ingraham, J.L., Low, K.B., Magasanik, B., Umbarger, H.E. (Eds.), Escherichia coli and Salmonella thyphymurium: Cellular and Molecular Biology, vol. 2. American Society for Microbiology, Washington, DC, pp. 1453–1472. Yanofsky, C., Horn, V., 1994. Role of regulatory features of the trp operon of Escherichia coli in mediating a response to a nutritional shift. J. Bacteriol. 176, 6245–6254. Zeron, E.S., 2008. Positive and negative feedback in engineering and biology. Math. Model. Nat. Phenom. 3 (2), 67–84.