Simulating in vitro transcriptional response of zinc homeostasis system in Escherichia coli.

BACKGROUND: The zinc homeostasis system in Escherichia coli is one of the most intensively studied prokaryotic zinc homeostasis systems. Its underlying regulatory machine consists of repression on zinc influx through ZnuABC by Zur (Zn2+ uptake regulator) and activation on zinc efflux via ZntA by ZntR (a zinc-responsive regulator). Although these transcriptional regulations seem to be well characterized, and there is an abundance of detailed in vitro experimental data available, as yet there is no mathematical model to help interpret these data. To our knowledge, the work described here is the first attempt to use a mathematical model to simulate these regulatory relations and to help explain the in vitro experimental data.
RESULTS: We develop a unified mathematical model consisting of 14 reactions to simulate the in vitro transcriptional response of the zinc homeostasis system in E. coli. Firstly, we simulate the in vitro Zur-DNA interaction by using two of these reactions, which are expressed as 4 ordinary differential equations (ODEs). By imposing the conservation restraints and solving the relevant steady state equations, we find that the simulated sigmoidal curve matches the corresponding experimental data. Secondly, by numerically solving the ODEs for simulating the Zur and ZntR run-off transcription experiments, and depicting the simulated concentrations of zntA and znuC transcripts as a function of free zinc concentration, we find that the simulated curves fit the corresponding in vitro experimental data. Moreover, we also perform simulations, after taking into consideration the competitive effects of ZntR with the zinc buffer, and depict the simulated concentration of zntA transcripts as a function of the total ZntR concentration, both in the presence and absence of Zn(II). The obtained simulation results are in general agreement with the corresponding experimental data.
CONCLUSION: Simulation results show that our model can quantitatively reproduce the results of several of the in vitro experiments conducted by Outten CE and her colleagues. Our model provides a detailed insight into the dynamics of the regulatory system and also provides a general framework for simulating in vitro metal-binding and transcription experiments and interpreting the relevant experimental data.

Background

Zinc is essential for life. It serves as a structural or catalytic cofactor in a large number of proteins such as RNA polymerase and zinc finger proteins [1-9]. Zinc also plays an important signalling role in various biological processes such as neurotransmission, cell proliferation, and apoptosis [10,11]. However, due to the potential toxicity of zinc, intracellular zinc concentrations must be kept under tight control. For example, a high intracellular Zn2+ concentration can inhibit the aerobic respiratory chain in E. coli [6-8].

E. coli achieves zinc homeostasis by regulating the uptake and efflux of zinc across the plasma membrane [1,8]. As we can see in Fig. 1a, extracellular zinc ions are transported into the cytoplasm through ZnuABC (an ABC-type transporter) and ZupT (a zinc permease), while the efflux of zinc is accomplished by ZntA (a P-type ATPase) and ZitB (a cation diffusion facilitator) [7,8,12-21]. Within the cytoplasm, similar to copper, it is thought that zinc trafficking may involve chaperone-like proteins [22,23]. However, despite considerable experimental effort, the zinc chaperone protein in E. coli has yet to be identified [2,6,24-27]. The ZnuABC transporter (encoded by the znuACB gene cluster) is composed of the periplasmic binding protein ZnuA, the ATPase ZnuC, and the integral membrane protein ZnuB [28]. This zinc uptake system is regulated by Zur, a dimer protein which binds at least 2 zinc ions. Zur is sensitive to the intracellular zinc concentration, and zinc-bound Zur (presumably the Zn4Zur form, the Zur dimer which contains 2 zinc ions per monomer and it is denoted as Zn2Zur in [22]) can compete with RNA polymerase to bind to the znu operator and act as a repressor [7,8,20].

Figure 1

Schematic representations of . (a) A schematic graph depicts the Zn2+ homeostasis system in Escherichia coli. Extracellular Zn2+ enters the cytoplasm through ZnuABC and ZupT [7,19]. In the presence of zinc, Zur binds to the znu operator and represses the transcription of znuACB gene cluster [8,20]. Excess intracellular zinc ions are exported by ZntA and ZitB [16,17,21]. Intracellular zinc can bind with protein ZntR and convert it into a strong transcriptional activator of the zntA gene [8,14,29]. The cytoplasmic zinc trafficking may involve chaperone-like proteins [22]. Abbreviations used in this graph are as follows: Zur* (active Zur); ZntR* (active ZntR); C? (zinc chaperone whose existence is still under debate) [2,22]. (b) A schematic graph depicts the main sub-processes which we need to model for simulating in vitro transcriptional response: (i) Zn2+-sensing by Zur, (ii) Transcriptional repression of znuC gene by Zur, (iii) Zn2+-sensing by ZntR and (iv) Transcriptional activation of zntA gene by ZntR (Please note that here we only model the transcription of znuC gene rather than of the whole znuACB gene cluster because we only have reported data for znuC transcripts available for comparison) [22].

In contrast to this mechanism, zinc efflux through ZntA is regulated by ZntR, a zinc-responsive MerR-like transcriptional regulator [8,14,29,30]. ZntR is a dimer protein which can bind one or two zinc ions per monomer depending on the buffer conditions [29]. A metal occupancy assay of ZntR, monitored by changes in tyrosine fluorescence, shows non-cooperative 1:1 binding of Zn(II) to the ZntR dimer [31]. ZntR in its apo form only slightly activates zntA transcription (please note that the apo form of ZntR (i.e., apo-ZntR) means that ZntR without the binding of Zn(II)) [8,14,29]. The binding of zinc-bound ZntR to the promoter introduces conformational changes in the DNA, which apparently make the promoter a better substrate for RNA polymerase, thus strongly activating the transcription of the zntA gene and increasing the efflux of zinc from the cell [29].

During 1999–2001, Outten CE and her colleagues presented some results on in vitro transcription and metal-binding competition experiments of E. coli zinc homeostasis system and showed that both ZntR and Zur are extremely avid zinc sensors and are both saturated at femtomolar free zinc concentrations [22,29,31]. In these experiments, the Zn(II) concentration was precisely controlled by using N,N,N',N'-tetrakis(2-pyridylmethyl) ethylenediamine (TPEN) as a zinc buffer [22]. The various assays relevant to this paper include the Zur-DNA interaction assay, Zur transcription assay and two ZntR transcription assays. In the Zur-DNA interaction assay the DNase I footprinting technique was used and the Zur-DNA interaction was found to correlate with the concentration of free Zn(II) (see the black dots in Fig. 2).

Figure 2

Simulation of Zur-DNA interaction. The black dots are reconstructed from the reported data in the original figure (the right graph in Fig. 3 in [22]) using image analysis method (please refer to the Methods for more details). The black curve is the simulated ratio (i.e., /D0 * 100%) of the final steady state concentration values of Zn4Zur-DNA complex (denoted by ) and the total concentration of DNA (D0 = 1 nM in this case) as a function of the logarithm of parameter Zn which denotes the simulated free zinc concentration.

In the Zur transcription assay, in vitro run-off transcription experiments with Zur and the znu Zn(II) uptake system were conducted, and the levels of the znuC RNA transcript were reported to correlate with the free Zn(II) concentration (see the red dots in Fig. 3d) [22]. In these real run-off transcription experiments, various reactants (including znuC DNA template, Zur, Zn(II), RNAP and heparin, etc.) were added sequentially and allowed to equilibrate first (~30 min total). Then nucleoside triphosphates (NTPs) were added and the reaction was stopped for 15 min (Outten CE, personal communication). Similar run-off transcription experiments (the ZntR transcription assay (I)) were conducted with ZntR and the zntA promoter and the levels of the zntA RNA transcript were reported to correlate with the free Zn(II) concentration (see the blue dots in Fig. 3d) [22]. Similarly in ZntR transcription assay (II), the levels of the zntA RNA transcript were reported to correlate with the total ZntR concentration, both with added Zn(II) and without Zn(II) (see the red and blue dots in Fig. 5a, respectively). Moreover, it was also reported that the levels of the zntA RNA transcript correlated with the total zinc concentration (see the black dots in Fig. 5b) [29].

Figure 3

Transient curves of simulated Zur transcription assay and data comparison (I). (a) The green, purple and red curves denote the simulated transient curves of Zn2Zur (Py), Zn4Zur (Py1), RNA polymerase (Rw) concentrations as a function of t, respectively. (b) The green, purple and red curves denote the simulated transient curves of free znuC DNA (Dw), znuC transcription initiation complex (Qw1) and Zn4Zur-DNA complex (Qw2) concentrations as a function of t, respectively. (c) The simulated concentration of the mRNA of ZnuC (Mw(t)) is depicted as a function of t. (d) Data comparison for Zur and ZntR transcription assays. Big red dots for the Zur transcription assay and big blue dots for the ZntR transcription assay (I) are reconstructed from the reported data in the original figure (Fig. 4 in [22]) using image analysis method. The purple curve and the green curve are the corresponding simulated normalized final concentrations of the mRNA of ZnuC (Mw (t = t+ t)) and the mRNA of ZntA (Mz (t = t+ t)) as a function of the logarithm of Zn (also denoted as [Zn]), respectively. The three small purple dots on the purple curve are simulated data points for Zn = 10-5 nM, 10-6 nM, 10-7 nM, respectively. The three small green dots on the green curve are simulated data points for Zn = 10-5 nM, 10-6 nM, 10-7 nM, respectively (please note that the simulated transient curves of Zur and ZntR transcription assays for Zn = 10-5 nM are shown in Fig. 3a-c and Fig. 4, respectively. More simulated transient curves for Zn = 10-6 nM, 10-7 nM are shown in Additional file 2: MoreTransientCurves.doc). The area highlighted in gray is the range of Zn between the half maximal induction points on the two simulated curves.

Figure 5

Comparison of simulated results and experimental data (II). (a) ZntR transcription assay with Zn(II) or without Zn(II). Red dots for the case of with Zn(II) and green dots for the case of without Zn(II) are reconstructed from the reported data in the original figure (Fig. 6B in [29]) using image analysis. Error bars indicate a standard deviation both above and below the average values of two separate experiments. The purple line and the green line are the corresponding simulated final concentrations (in nM) of mRNA of ZntA (Mz(t = t+ t)) in the cases of parameter Zn= 10 μM and Zn= 0 as a function of the logarithm of [ZntR](i.e., 2*Px), respectively. (b) ZntR transcription assay with varying total zinc concentration. The black dots are reconstructed from the reported data in the original figure (Fig. 6C in [29]) using image analysis. The black curve is the simulated final concentration (in nM) of mRNA of ZntA (Mz(t = t+ t)) as a function of the logarithm of parameter Zn(also denoted as [Zn]).

Although the transcriptional regulation of the zinc homeostasis system in E. coli seems to be well characterized, and despite the fact that detailed in vitro experimental data on this system are also available, as yet there is no mathematical model to help interpret these data [22,29,31]. The principal aim of this paper is to present a mathematical model which is capable of simulating this regulatory system and can be used to help interpret various experimental data.

We will present a unified mathematical model and use it to simulate the in vitro transcriptional response of the zinc homeostasis system in E. coli. The construction of the model is based on biochemical principles and we use open source software (Cellerator) to automatically generate the equations [32,33]. We validate our model by comparing the simulation results with the corresponding in vitro experimental data.

Results

As shown in Table 1, we use 14 reactions to represent the four sub-processes involved in the zinc homeostasis system, namely: (i) Zn2+-sensing by Zur, (ii) transcriptional repression of the znuC gene by Zur, (iii) Zn2+-sensing by ZntR and (iv) transcriptional activation of the zntA gene by ZntR (see Fig. 1b, and for more details please refer to the Methods). Here we present our results for simulating various in vitro assays (please note that the main differences between the ZntR transcription assay (I) and assay (II) are differences in the initial conditions and in that in assay (II) we take into consideration competition between ZntR and TPEN for zinc binding by including Reaction 14, whereas in assay (I), Reaction 14 is not included).

Table 1

The reactions of the model

Sub-Process Name	Reaction No.	Cellerator Form of Particular Reactions	Description
Zn2+-Sensing by ZntR	(1)	{Px + Zn ⇄ Px1, r1, r2}	apo-ZntR binding with zinc to become active ZntR
Transcriptional Activation of zntA Gene by	(2)	{Dz + Rz ⇄ Qz1, k2a, k-2}	DNA of ZntA binding with RNAP
	(3)	{Qz1 ⇄ Dz + Mz + Rz, k3, 0}	transcription of Qz1
	(4)	{Dz + Px ⇄ Qz4, k1b, k-1}	apo-ZntR binding with DNA
	(5)	{Qz4+ Rz ⇄ Qz5, k2b, k-2}	apo-ZntR-DNA complex binding with RNAP
	(6)	{Qz5 ⇄ Qz4+ Mz + Rz, k3, 0}	transcription of Qz5
	(7)	{Dz + Px1 ⇄ Qz2, k1, k-1}	ZnZntR binding with DNA
	(8)	{Qz2 + Rz ⇄ Qz3, k2c, k-2}	ZnZntR-DNA complex Binding with RNAP
	(9)	{Qz3 ⇄ Qz2+ Mz + Rz, k3, 0}	transcription of Qz3
Zn2+-Sensing by Zur	(10)	{Zn2+ Py ⇄ Py1, r3, r4}	Zn2Zur binding with zinc to become active Zur
Transcriptional Repression of znuC Gene by Zur	(11)	{Dw + Py1 ⇄ Qw2, k1a, k-1}	active Zur binding with DNA to form complex Qw2 which can not bind with RNAP
	(12)	{Dw + Rw ⇄ Qw1, k2, k-2}	DNA of ZnuC binding with RNAP
	(13)	{Qw1 ⇄ Dw + Mw + Rw, k3, 0}	transcription of Qw1
Zn2+-Binding by TPEN	(14)	{Zn + Tp ⇄ Tp1, r5, r6}	TPEN binding with zinc to form a complex

Note: Abbreviations and synonyms used in this table are as follows: Zn (free zinc ion); Px (apo-ZntR); Px1 (active ZntR, i.e., ZnZntR); Py (the Zur dimer which contains two zinc ions per dimer, here we denote it as Zn2Zur and it is denoted as Zn1Zur in [22]); Py1 (active Zur, i.e., the Zur dimer which contains four zinc ions per dimer, here we denote it as Zn4Zur and it is denoted as Zn2Zur in [22]); Z (ZntA); DZ (DNA of ZntA); Rz (RNA polymerase for zntA transcription); Mz (mRNA of ZntA); Qz1 (transcription initiation complex formed by Dzand Rz); Qz2 (ZnZntR-DNA complex); Qz3 (transcription initiation complex formed by Qz2 and Rz); Qz4 (apo-ZntR-DNA complex); Qz5 (transcription initiation complex formed by Qz4 and Rz); w (ZnuC);Dw (DNA of ZnuC); Rw (RNA polymerase for znuC transcription); Mw (mRNA of ZnuC); Qw1 (transcription initiation complex of ZnuC); Qw2 (Zn4Zur-DNA complex which can not further bind with Rw); Tp (free TPEN not bound by zinc); Tp1 (zinc-bound TPEN); RNAP (RNA polymerase).

Zur-DNA interaction

The Zur-DNA interaction assay involves only two reactions (Reactions 10 and 11, see Table 1), which are expressed as 4 ODEs (for the detailed equations, see Additional file 1: ModelEquations.doc) [22]. By imposing the conservation restraints (Py(t) + Py1(t) + Qw2(t) = Py= 25 nM, Dw(t) + Qw2(t) = D0 = 1 nM) (as in the real experiment [22]) and solving the relevant steady state equations (for the parameters, please see Table 2. Note that in the real experiment the total concentration of Zur monomer is 50 nM, here we need to divide this value by half which means that Py= 25 nM because in solution, Zur exists in dimer form [20]. In similar way we can calculate Px), we can depict the simulated ratio of steady state concentrations of the Zn4Zur-DNA complex (denoted by ) and the total concentration of znuC DNA (D0 = 1 nM in this case) as a function of the logarithm of parameter Zn as shown in Fig. 2 (the black curve). From this figure, we can see that when the simulated free zinc concentration (Zn) ranges from 10-18 M to 10-14 M, the simulated protection ratio (denoted by /D0 * 100%) rises from 0.00014% to 96.4%. This means that in the presence of higher free zinc concentrations, more Zn2Zur molecules become active and bind with znuC DNA molecules to protect them from the binding of RNA polymerase. The simulated sigmoidal curve (the black curve in Fig. 2) seems to fit well with the corresponding experimental data (the black dots in Fig. 2) [22].

Table 2

Model parameters for which all results are calculated unless otherwise stated

Parameter	Value	Description
k d	10-14.9 M	the Zn(II) dissociation constant for ZnZntR when pH = 8.0 [31]
kd1	10-15.2 M	the Zn(II) dissociation constant for the ZnZntR-DNA complex when pH = 8.0 [31]
K Z n − T P E N '	1.99*1015 M-1	the apparent association constant for Zn-TPEN at pH = 8.0, 0.1 M ionic strength, calculated from [22]
k1	0.025 (nM)-1s-1	the forward rate parameter of Reaction (7)
k1a	1 (nM)-1s-1	the forward rate parameter of Reaction (11)
k1b	1.253*10-2 (nM)-1s-1	the forward rate parameter of Reaction (4)
k-1	0.9 s-1	the backward rate parameter of Reactions (4,7,11)
k2	0.02 (nM)-1s-1	the forward rate parameter of Reaction (12)
k2a	0.00005 (nM)-1s-1	the forward rate parameter of Reaction (2)
k2b	0.0002 (nM)-1s-1	the forward rate parameter of Reaction (5)
k2c	0.0037 (nM)-1s-1	the forward rate parameter of Reaction (8)
k-2	0.3 s-1	the backward rate constant of Reactions (2, 5, 8, 12)
k3	0.011 s-1	the transcription rate parameter
r1	2.73*102 (nM)-1s-1	the forward rate parameter of Reaction (1)
r2	3.437*10-4 s-1	the backward rate parameter of Reaction (1)
r3	4.41*1010 (nM)-2s-1	the forward rate parameter of Reaction (10)
r4	9*10-3 s-1	the backward rate parameter of Reaction (10)
r5	3*104 (nM)-1s-1	the forward rate parameter of Reaction (14)
r6	1.506 *10-2 s-1	the backward rate parameter of Reaction (14)
td0	30 min	the time duration for preliminary equilibrium of reactants before NTPs (i.e., nucleoside triphosphates)were added in run-off transcription experiments [22,29,37]
t d	15 min	the time duration for run-off transcription after NTPs were added in transcription experiments [22,29,37]
td1	30 min	the time duration for Zur-DNA interaction assay [22]
Px tot	25 nM	the total concentration of ZntR dimer which is half of the concentration of ZntR monomer denoted as [ZntR]total[22]
Py tot	25 nM	the total concentration of Zur dimer [22]
R0	50 nM	the total concentration of RNAP [22]
D0	4 nM	the total concentration of DNA [22]
Zn tot	vary in different assays	the total concentration of Zn(II), also denoted as [Zn]total
TPEN tot	vary in different assays	the total concentration of TPEN

Note: k2= 4*k2. Moreover, according to the equilibrium theory of chemical reactions, r2 = kr1, r6 = r5/ and k1= k1k/k. The values of four parameters (k1, k-1, k2, k-2) are taken from Hayot's model [43]. These parameters origin from measured rate constants of the λ repressor gene cI in E. coli and are also quoted as physiologically reasonable values by Ingram et al [15,43,44]. is calculated in the same way as shown in [22] (please note the pH value difference).

We derive the same simulation results by directly solving the 4 relevant ODEs with Py(0) = Py= 25 nM, Dw(0) = D0 = 1 nM, Dw2(0) = 0, Py1 (0) = 0 as the initial conditions and depicting the simulated ratio of the final concentration of the Zn4Zur-DNA complex (Qw2(t = t) and D0. This is because the system reaches equilibrium before t = t= 30 min.

Zur transcription assay

As mentioned in the legend of Fig. 1b, here we only simulate the transcription of the znuC gene. We approximate the in vitro Zur run-off transcription assay by a two-phase (namely, the preliminary equilibrium phase and the transcription phase) sub-model. In the first phase, the preliminary equilibrating process of reactants involves 3 reactions (Reactions 10–12) which are expressed as 6 ODEs (see Additional file 1: ModelEquations.doc). In the second phase, the run-off transcription involves 4 reactions (Reactions 10–13 because now the real transcription happens after the addition of the NTPs) which are expressed as 7 ODEs (see Additional file 1: ModelEquations.doc).

By setting the initial conditions of the model simulation to be the same as those in the real experiment (Py(0) = Py= 25 nM, Dw(0) = D0 = 4 nM, Rw(0) = R0 = 50 nM and all the remaining initial concentrations are set to be 0) and numerically solving the 6 equations for the first phase and then solving the 7 ODEs for the second phase (obviously we need to use the end concentration values of the reactants in the first phase as the initial concentration values of reactants in the second phase), we can depict the relevant transient curves for Zn = 10-5 as shown in Fig. 3a–c (for the values of the remaining parameters, please see Table 2).

As shown in Fig. 3a, due to the binding of free zinc, the simulated concentration of Zn2Zur (Py(t), the Zur dimer which contains two zinc ions per dimer and it is used in the corresponding real assay [22]) quickly decreases from 25 nM to a steady state value of 0.044 nM whereas the simulated concentration of active Zur (Py1 (t), the Zur dimer which contains four zinc ions per dimer) quickly rises from 0 to 21.6 nM. The simulated concentration of RNA polymerase (Rw(t)) decreases slightly from 50 nM to 49.5 nM due to the effect of its binding with znuC DNA.

As we can see from Fig. 3b, the simulated free znuC DNA concentration (Dw(t)) decreases rapidly (in 0.4 seconds) from 4 nM to a steady state of 0.14 nM during the first phase due to the binding of active Zur and RNA polymerase. The simulated concentration of the transcription initiation complex (Qw1(t)) rapidly increases (in 0.4 seconds) from 0 to a steady state value of 0.47 nM whereas the simulated concentration of Zn4Zur-DNA complex (Qw2(t)) quickly increases (in 0.6 seconds) from 0 to 3.39 nM. The initiation of the second phase seems to only have a small influence on the afore mentioned steady state values (e.g., the steady state values of Qw1(t) and Qw2(t) change from 0.47 nM and 3.39 nM at the end of first phase to 0.45 nM and 3.4 nM at the end of the second phase, respectively). From Fig. 3c, we can see that in the first 30 minutes, the concentration of mRNA of ZnuC (Mw(t)) remains at 0 because the real transcription has not happened yet, and then in the subsequent 15 minutes it increases linearly from 0 to a final concentration of 4.49 nM.

The rapid decrease in the concentration of free znuC DNA (Dw(t)) shown in Fig. 3b is due to the binding of znuC DNA with active Zur (Zn4Zur) and RNA polymerase. Since in the whole process, the total increase in the simulated concentration of Zn4Zur-DNA complex (Qw2(t)) is 3.4 nM, whereas the total decrease of the simulated free znuC DNA concentration is about 3.86 nM, we can conclude that when Zn = 10-5 nM, the binding of active Zur consumes the majority of the znuC DNA to form the Zn4Zur-DNA complex, which can not further bind with RNA polymerase, and in this way the transcription of znuC is repressed.

We performed many simulations for various values of Zn (in the range of 10-18 M to 10-12 M) and recorded the final values of the simulated mRNA concentration (Mw(t = t+ t)). After normalizing these concentration values, depicting them as a function of Zn (in logarithm), and smoothly connecting these simulated data points, we obtained the purple curve in Fig. 3d (please note that only three simulated data points for Zur assay are shown as small purple dots in this figure to avoid confusion with the experimental data points).

ZntR transcription assay (I)

Similarly, the ZntR run-off transcription assay can also be simulated by a two-phase sub-model. The first phase (the preliminary equilibrium phase) involves 6 reactions (Reactions 1,2,4,5,7,8), which are expressed as 9 ODEs (see Additional file 1: ModelEquations.doc). The second phase (the transcription phase) involves 9 reactions (Reactions 1–9), which are expressed as 10 ODEs (see Additional file 1: ModelEquations.doc). By setting the initial conditions of the simulation to be the same as those used in the real experiment (Px(0) = Px= 25 nM, Dz(0) = D0 =4 nM, Rz(0) = R0 = 50 nM and setting all the remaining initial concentrations to be 0), and subsequently solving the relevant equations of the two-phase sub-model, we can depict the relevant transient curves for Zn = 10-5 nM as shown in Fig. 4 (for the remaining parameters, please see Table 2). In this assay, Reaction 14 is not included in the sub-model because the ZntR concentration is too low to challenge the buffering capacity of TPEN (of course we can also perform numerical simulations by including Reaction 14, although further investigations have shown that we essentially get the same results).

Figure 4

Transient curves of simulated ZntR transcription assay (I). (a) The green, purple and red curves denote the simulated transient curves of apo-ZntR (Px), ZnZntR (Px1), RNA polymerase (Rz) concentrations as a function of t, respectively. (b) The simulated concentrations of the free zntA DNA (Dz, green curve) and mRNA of ZntA (Mz, purple curve) are depicted as a function of t. (c) The green and purple curves denote the simulated transient curves of transcription initiation complexes (Qz1 and Qz5) concentrations as a function of t, respectively. (d) The green, purple and red curves denote the simulated transient curves of ZnZntR-DNA complex (Qz2), transcription initiation complex (Qz3) and apo-ZntR-DNA complex (Qz4) as a function of t, respectively.

From Fig. 4a we can see that due to the binding with free zinc, the simulated concentration of apo-ZntR (Px(t)) decreases from 25 nM to 2.59 nM, whereas the simulated concentration of active ZntR (Px1(t)) rises from 0 to 20.5 nM, and the simulated concentration of RNA polymerase (Rz(t)) decreases slightly from 50 nM to 49.3 nM. In the first phase, due to the binding with ZntR and RNA polymerase, the simulated unbound zntA DNA concentration (Dz(t)) decreases rapidly (in 0.04 minutes) from 4 nM to 3.03 nM and then decreases gradually to 2.04 nM at the end of the first phase (Fig. 4b, green curve); in the second phase, the free zntA DNA concentration remains at roughly the same level (2.05 nM). The simulated zntA mRNA concentration (Mz(t)) remains at 0 nM during the first phase, as there is no transcription happening, and then increases seemingly linearly to a final concentration of 6.96 nM during the second phase after NTPs have been added (Fig. 4b, purple curve).

The simulated transients curves in Fig. 4c show that Qz1(t) rapidly rises (in 0.3 minutes) from 0 to a peak value of 0.024 nM and then gradually decreases to 0.017 nM during the first 30 minutes whereas Qz5(t) rapidly rises (in 0.26 minutes) from 0 to a peak value of 0.031 nM and then gradually decreases to a 0.0025 nM during the first phase. The initiation of the second phase causes a small decrease in the values of Qz1(t) and Qz5(t) (to 0.016 nM and to 0.0023 nM, respectively). As shown in Fig. 4d, both Qz2(t) and Qz3(t) rise first (from 0 to 1.16 nM and 0.7 nM, respectively) during the first phase whereas Qz4(t) first dramatically increases up to a peak value of 0.97 nM and then gradually decreases to its final value of 0.076 nM. The initiation of the second phase causes a small decrease in the value of Qz3(t) and a slight increase in the value of Qz2(t), as judged by the small kinks in the corresponding two curves, whereas it has insignificant influence of the value of Qz4(t).

Using similar methods we can obtain the green curve in Fig. 3d for the final values of the simulated zntA mRNA concentration (Mz(t = t+ t)) as a function of the value of Zn (in logarithm). The results shown in Fig. 3d indicate that when the simulated free zinc concentration ranges from 10-18 M to 10-12 M, the simulated normalized final concentrations of mRNA of ZnuC (Mw(t = t+ t)) decreases from 100% to 15.05%, whereas the simulated normalized final concentration of mRNA of ZntA (Mz(t = t+ t)) increases from 7.4% to 100%. The half-maximal induction of znuC transcripts and the half maximal induction of zntA transcripts occur at Zn = 2*10-16 M and Zn = 1.15*10-15 M respectively, which are the same as previously reported values [22]. The simulated purple curve (for the Zur transcription assay) agrees with corresponding experimental data (the red dots) extremely well. Similarly the simulated green curve (for the ZntR transcription assay (I)) also agrees with the corresponding experimental data (the blue dots), although to a slightly lesser degree [22].

ZntR transcription assay (II)

In this assay, we take into consideration the competition between ZntR and TPEN for zinc binding by including Reaction 14. Again, we will use a two-phase sub-model to simulate the real assay. The first phase (the preliminary equilibrium phase) of the assay (II) involves 7 reactions (Reactions 1,2,4,5,7,8,14), which are expressed as 12 ODEs (see Additional file 1: ModelEquations.doc). The second phase (the transcription phase) involves 10 reactions (Reactions 1–9,14), which are expressed as 13 ODEs (see Additional file 1: ModelEquations.doc). By setting the initial conditions of the simulation equal to those used in the real experiment (, and all the remaining initial concentrations are set to 0) and solving the two-phase model, we depict the simulated final concentrations (in nM) of mRNA of ZntA (Mz(t = t+ t)) for Zn= 10 μM and Zn= 0 as a function of the logarithm of the doubled value of parameter Px(i.e., [ZntR]which denotes the total concentration of the ZntR monomer) and we obtain the purple and green curves shown in Fig. 5a[29].

We also perform many simulations under the following initial conditions: (, and all the remaining initial concentrations are set to 0) for various values of Znwithin the range of 100 nM to 100 μM and eventually obtain the black curve shown in Fig. 5b which describes the final values of the simulated mRNA concentration (Mz(t = t+ t)) as a function of the value of Zn(in logarithm).

Discussion

The simulation results shown in Fig. 4 indicate the complex interactions among three transcription processes of zntA (the constitutive transcription, the apo-ZntR activated transcription and the ZnZntR activated transcription). If we compare the dynamics of the simulated concentrations of three transcription initiation complexes involved in the ZntR transcription assay (i.e., Qz1(t), Qz3(t) and Qz5(t)) as shown in Fig. 4c and Fig. 4d, we find that the dynamics of Qz1(t) and Qz5(t) are quite similar. Initially, they both increase rapidly, form low peaks (the peak values are 0.024 nM and 0.031 nM, respectively), and then gradually decrease. In contrast, the dynamics of Qz3(t) only demonstrates a gradual increase to 0.7 nM in the first 30 minutes. The observation that the final steady state value of Qz3(t) (0.69 nM) is much higher than those of Qz1(t) and Qz5(t) (0.016 nM and 0.0023 nM, respectively) indicates that for Zn = 10-5nM, when the system (excluding Mz(t)) enters its final equilibrium, the dominating transcription process is ZnZntR activated transcription rather than the other two transcription processes (i.e., the constitutive transcription and the apo-ZntR activated transcription, please refer to Methods for more details).

To explain why the dynamics of Qz1(t) shows a peak, we suggest that the initial increase of Qz1(t) is due to the binding of zntA DNA with RNA polymerase. Then following the conversion of apo-ZntR to active ZntR by zinc-binding (see the green and purple curves in Fig. 4a), active ZntR binds with zntA DNA to form the ZnZntR-DNA complex (see the green curve in Fig. 4b and the green curve in Fig. 4d). This competitive binding of active ZntR causes a sudden decrease in the free zntA DNA concentration (see the green curve in Fig. 4b) and the reversible Reaction 2 (see Table 1) becomes dominated by its reverse side and Qz1(t) begins to decrease after forming a small peak. Similarly, we can explain the dynamics of Qz5(t).

By comparing the dynamics of the simulated Zur and ZntR transcription assays shown in Fig. 3a–c and Fig. 4, we can see that when Zn = 10-5nM, during the first phase, the simulated Zur transcription system reaches its steady state in less than 20 seconds, much faster than the simulated ZntR transcription system which takes more than 20 minutes. As shown in Fig. 3c and Fig. 4b, the seemingly linear increase of the simulated concentrations of mRNA (Mw(t) and Mz(t)) during the second phase indicates the progress of the relevant transcription processes. If we calculate the slope of the linear curve in Fig. 3c as follows:

Mw(t = t+ t)//t= 4.49 nM/0.45 nM/15 min = 0.011s-1 where denotes the final steady state value of Qz1(t), we derive the same value as that of the transcription rate parameter k3. Obviously the simulated final concentrations of mRNA (Mw(t = t+ t) and Mz(t = t+ t)) are generally proportional to t, which is in accordance with the experimental observation that the harvest of run-off transcription assay is related to the duration time of its transcription phase (t) [34].

The purple curve in Fig. 5a indicates that for Zn= 10 μM, when the simulated total ZntR monomer concentration ([ZntR]which is twice the value of Px) ranges from 10-9 M to 10-4 M, the simulated final concentrations of mRNA of ZntA (Mz(t = t+ t)) increases from 0.59 nM to a peak value of 10.09 nM when [ZntR]= 10-5.92M and then decreases to 4.99 nM. If we look at the corresponding experimental data (the red dots), we can see that the relative induction of the zntA transcripts increases, forms a peak (when [ZntR]= 10-7M), and eventually declines [29]. Thus our simulation successfully simulates the peak behaviour of the relative induction of the zntA transcripts for increasing values of [ZntR]in the presence of zinc. Further investigations show that if we perform the simulations excluding Reaction 14, then we can only reproduce the increasing behaviour rather than the peak behaviour. Thus one potential explanation for the peak behaviour is that, for low ZntR concentrations, TPEN is strong enough to buffer the zinc and more ZntR will promote the transcription of zntA gene; while for high ZntR concentrations, the buffering capacity of TPEN is exceeded and the free zinc concentration can not be maintained as a constant anymore and it subsequently decreases due to the binding of over-abundant ZntR molecules, which in turn limits the transcription processes. A similar comparison can be made for the case when Zn= 0 (i.e., in the absence of zinc, please see the green curve and the blue dots in Fig. 5a). However, in the latter case, our model can only simulate the initial increase, but fails to reproduce the decline.

As described in detail in the Methods section, in this model, we assume that the active form of ZntR is ZnZntR because metal occupancy assays of ZntR monitored by changes in tyrosine fluorescence show noncooperative 1:1 binding of Zn(II) to the ZntR dimer [31]. This assumption is valid only when the free zinc concentration and total ZntR concentration are both extremely low (in sub-nM and nM range, respectively). When the total ZntR concentration goes to the μM range, the binding kinetics of Zn(II) to the ZntR dimer will be more complicated because ZntR can bind one or two zinc ions per dimer depending on the buffer conditions [8,29]. This explains why, as shown in Fig. 5a in the case of with Zn(II), there is a disagreement between the simulation results (the purple curve) and the corresponding experimental data (the red dots) when ZntR molecules are relatively abundant. Intuitively, we can think of it in this way: in the real case, the competitive ability of ZntR for Zn(II) binding is stronger than the model prediction because at high ZntR concentrations, ZntR, on average, binds with more than one ion per dimer. This results in a smaller and earlier peak because the buffering capacity of TPEN is now easier to exceed. In the absence of Zn(II), the eventual abnormal decline in the experimental data (see the blue dots in Fig. 5a) may be due to the normal deviations of the different experiments because the levels of zntA transcript are very low in this case or perhaps this is due to some novel, unknown mechanisms (please note that the error bars shown in Fig. 5a indicate the standard deviation from the average values of only two separate experiments and there are only two data points having error bars for the case of without Zn(II)) [29].

As we can see from Fig. 5b, when the simulated total zinc concentration (Zn) ranges from 10-7 M to 10-4 M, the simulated final concentration of mRNA of ZntA (Mz(t = t+ t)) increases from 0.72 nM to 8.04 nM (saturation occurs when Zn= 10-5M) which means that more abundant free zinc ions bind with ZntR to promote activation of the transcription of the zntA gene. The simulated curve (the black curve) fits the experimental data (the black dots) quite well [29].

Coupled feedback loops have been recently recognized as essential building blocks of cellular networks [35]. The zinc homeostasis system in E. coli is a good example of such a building block because it follows from Fig. 1a that Zur and ZnuC form a 'negative circuit', since active Zur represses znuC (negative action) while zinc influx via ZnuC leads to larger amounts of active Zur molecules (positive action). Similar considerations point towards the negative circuit wiring between ZntR and ZntA. It is believed that such coupled negative feedback loops are quite helpful for enhancing homeostasis [35].

As previously mentioned, cytoplasmic zinc trafficking in E. coli may involve chaperone-like proteins whose existence is still being debated [2,6,22]. Outten et al. demonstrated in vitro that ZntR and Zur are sensitive to very low concentrations (femtomolar) of free zinc (also see Fig. 3d), therefore they proposed that free zinc in the cytosol of E. coli is not physiologically available under normal growth conditions [22]. Our simulation results further confirm their experimental data and support their proposal. However, in order to better understand the in vivo transcriptional regulation mechanisms of zinc homeostasis, further investigations are required to simulate the in vivo transcription processes and their responses to various environmental conditions.

Up until now, performing well-designed in vitro experiments has been one of the common ways used to infer the various characteristics of the corresponding in vivo systems. The current work provides a good example of how to use a unified mathematical model to explain complicated datasets obtained from in vitro metal-binding and transcription experiments which have been widely performed for metal ion homeostasis and detoxification systems [22,29,31,36,37]. The repression of Zur on the transcription of znuACB gene cluster and the activation of ZntR on the zntA transcription constitute the critical parts of the regulatory mechanisms of the zinc homeostasis system in E. coli (see Fig. 1a). This means that if we want to make predictive and useful model for the in vivo zinc homeostasis system, we need to model these transcriptional regulations. Although the current model only simulates the in vitro kinetics, together with its fitted rate constants it can be used as a good basis and reference for the future modelling of the corresponding in vivo system. Moreover, the quantitative distinguishment of the three transcription processes of zntA (the constitutive transcription, the apo-ZntR activated transcription and the ZnZntR activated transcription) in our model will be quite meaningful for modelling the in vivo system and it provides the possibility of including any additional regulations on these three processes which do happen in vivo [14,29].

In order to further our understanding of the process of zinc homeostasis in E. coli, the most critical thing is to identify the intracellular zinc chaperone, which is very likely to exist. Recently, proteomics has progressed to such a stage that it can determine the cellular response to any perturbation at the level of protein activation [38-40]. Thus mass spectrometry-based proteomics can be used to search the possible molecular candidates in addition to genome-wide high-throughput screens [25]. Once the zinc chaperone has been identified, the next step will be to measure the interactions between the zinc chaperone and the membrane transport proteins (ZnuABC, ZupT, ZntA and ZitB) and the interactions between the zinc chaperone and the metalloregulatory proteins (ZntR and Zur). Since similar work has already been done for the copper homeostasis system in E. hirae, the same equipments and experimental techniques used there (e.g., surface plasmon resonance analysis) can also be used to measure the kinetics of these interactions in the zinc homeostasis system [23,41].

Furthermore, we need to further take into consideration zinc storage and zinc utilisation by proteins in E. coli and quantify the concentrations of the relevant proteins, DNAs and mRNAs. Finally, the subtle details of relevant regulatory processes (e.g., proteolysis which has been proven to play a role) need to be further investigated [8]. Once we have characterised these processes in detail and have made corresponding sub-models for them, we can then integrate these sub-models together with the current model, in order to build a comprehensive model to describe the entire in vivo system. Further experiments determining the in vivo cellular response to various perturbations will be necessary for checking the validity of the model and also for model refinement. In this way, step by step we will acquire a complete map of the zinc homeostasis system in E. coli and reach a full understanding of the system dynamics. Close cooperation between pioneering experimentalists and computational scientists through iterative systems biology procedure (model → experiments → model) will be necessary for achieving such ambitious goals [42].

Conclusion

To summarize, we have built a mathematical model for simulating the in vitro transcriptional response of zinc homeostasis system in E. coli. Simulation results show that our model can quantitatively reproduce the various results of the in vitro experiments conducted by Outten CE and her colleagues. Our model gives a detailed insight into the involved system dynamics and provides a general framework for simulating in vitro metal-binding and transcription experiments and interpreting relevant experimental data.

Methods

Cellerator software

Cellerator™ is a Mathematica® package designed to facilitate biological modeling via automated equation generation [32,33]. It uses an arrow-based reaction notation to represent biochemical networks and is especially amenable for simulating signal transduction networks. For example, a reversible biochemical reaction (A + B ⇌ C, which means reactant A binds with reactant B to form product C, can be represented as {A + B ⇄ C, r, r} in Cellerator form where rand rdenote the forward and the backward rate constants, respectively. The detailed ODE notation of this reaction is: .

Representation of relevant reactions

As we can see in Fig. 1b, in order to simulate the in vitro transcriptional response, we need to model the four sub-processes involved:

1) Zn2+-sensing by ZntR

ZntR is a dimer protein which can bind one or two zinc ions per dimer depending on the buffer conditions [8,29]. However, metal occupancy assay of ZntR monitored by changes in tyrosine fluorescence shows non-cooperative 1:1 binding of Zn(II) to the ZntR dimer [31]. The zinc-bound form of ZntR has been reported to contain 0.75 ± 0.075 zinc/monomer, neither favoring 1:1 binding nor 1:2 binding [8]. However, this result was obtained under the condition of excessive ZntR protein (5 μM) [8]. Since the free zinc concentration and total ZntR concentration are both extremely low (in sub-nM and nM range, respectively) in all the relevant real assays (except the ZntR transcription assay (II) related to Fig. 5a) of this paper, here we assume that the active form of ZntR is ZnZntR (i.e. there is a 1:1 binding) [22,29]. We use Reaction (1) (see Table 1) to describe this sub-process.

2) Transcriptional activation of zntA gene by ZntR

Experimental results have shown that there is constitutive transcription activity of the zntA promoter [14]. According to Hayot et al., this constitutive transcription can be described by Reactions (2–3) (see Table 1, please note that the justification for the specific parameter values used in Hayot's model can be found in [15]. Hayot's model is later used by Ingram et al. to study the dynamics of the bi-fan motif) [43,44].

In the absence of Zn(II), apo-ZntR binds to the promoter and distorts the DNA which appears to result in an approximately fourfold induction [14]. According to Hayot et al., this apo-ZntR activated transcription can be described by Reactions (4–6) (see Table 1) and we have the relation: k2= 4* k2[43].

The binding of Zn(II) to ZntR converts it into a transcriptional activator protein that introduces conformational changes in the DNA which apparently make the promoter a better substrate for RNA polymerase [29]. According to Hayot et al., this ZnZntR activated transcription can be described by Reactions (7–9) (see Table 1) [43].

3) Zn2+-sensing by Zur

Zur is a dimer protein which binds at least 2 zinc ions [20,22]. Experimental results have established that the DNA binding of Zur presumably involves the Zn4Zur form (i.e., the Zur dimer which contains 2 zinc ions per monomer and it is denoted as Zn2Zur in [22]) rather than the Zn2Zur form (the Zur dimer which contains one zinc ion per monomer and it is denoted as Zn1Zur in [22]). Similar as Cui et al. did for modelling the binding of calmodulin with calcium ions, we use Reaction (10) (see Table 1) to describe this sub-process under the assumption of strong cooperativity existing between the two active sites of Zn2Zur (please note that the purified Zur dimer which contains one zinc ion per monomer is used in the relevant assays) [22,45].

4) Transcriptional repression of znuC gene by Zur

The genes znuA and znuCB are transcribed divergently and both promoters of znuA and znuCB are active in vivo [7,22]. Since we only have reported data for znuC transcripts available for comparison, here we choose to model the transcription of the znuC gene only [22]. In the absence of Zn(II), Zur does not compete for DNA binding. The addition of excessive Zn(II) allows Zur to bind to the znuC promoter and prevents its binding with RNA polymerase [22]. According to Hayot et al., we can use Reactions (11–13) (see Table 1) to describe this process [43].

Zinc binding by TPEN

As mentioned before, TPEN is used as a zinc buffer to precisely control the free zinc concentration in the relevant assays and this process can be apparently described by Reaction 14 [22,29]. Normally the free zinc concentration (Zn) is regarded as a constant and it can be simply calculated from the total zinc concentration (Zn) according to the following buffer equation:

However, in more complicated cases such as the ZntR transcription assay (II), it is wiser to perform numerical simulations by including this reaction and the free zinc concentration is no longer regarded as a constant.

The equations of the model and the numerical solver

The detailed equations used for simulating different assays can be found in Additional file 1: ModelEquations.doc. We use Mathematica's differential equation solver "NDSolve" to solve the relevant ODEs. If the studied ODEs are stiff as is the case for the relevant simulations of Fig. 5, we set the method option of NDSolve to be "StiffnessSwitching".

Translating the model into CellML

CellML is an XML-based modelling language which provides an unambiguous method of defining models of biological processes [46,47]. The current model has been translated into two CellML versions [48,49]. The first version (please visit the webpage for downloading the detailed code) is for ZntR transcription assay (I) which excludes the buffering equation of TPEN (i.e., Reaction 14) [48]. The second version (please visit the webpage for downloading the detailed code) is for ZntR transcription assay (II) which includes the buffering reaction of TPEN [49].

The image analysis method

The original figures are imported into the Paint tool of Windows system. The pixel coordinates are recorded for the axis origin, two tick points (one tick point on the horizontal axis and one tick point on the vertical axis) and all experimental data points. Then by simple algebraic calculations we can get the real coordinate values of the reported data points. For example, imagine that we need to analyze an image with x coordinate (in logarithm) and normal y coordinate. Assume the measured pixel coordinates of the axis origin (its real coordinate values are {10a, b}) and tick points (their real coordinate values are {10a, c} and {10d, b}) are (px0, py0), (px0, py1), (px1, py0), respectively. For a data point with measured pixel coordinates (px2, py2), we can calculate its real coordinate values {f, g} as follows:

The relative error of such data reconstruction is estimated to be (0.5–3)% depending on the image size.

Authors' contributions

JC conceived of the study, participated in its design, performed numerical simulations and helped to draft the manuscript. JAK participated in its design and helped to draft the manuscript. CML translated the model into CellML versions and helped to draft the manuscript. All authors read and approved the final manuscript.

Additional file 1

ModelEquations.doc. This additional file describes the detailed equations for simulating various assays.

Click here for file

Additional file 2

MoreTransientCurves.doc. This additional file describes the simulated transient curves of Zur and ZntR transcription assays for parameter Zn = 10-6nM, 10-7nM.

Click here for file