Modelling the series of reactions relating Iron (II) ion concentration in saliva to amilCP concentrations.

Abstract

This year, our team created a mathematical model to represent the series of reactions from Iron(II) ion concentrations in saliva to amilCP concentrations for our team’s system. We first had to decide whether we should model our final cell-free system (in vitro), or our in vivo lab work. An in vivo model would be advantageous, as it would allow us to compare the lab work we did to to a model. However, we decided an in vitro model would be more appropriate. This is because the in vivo lab work was just to test whether our construct worked, not to find a relationship between iron (II) ions and amilCP. We therefore felt a model to represent the relationship in an in vitro environment would better represent the final product our project aimed to create.

It should be noted that in a cell-free system, we would have to provide the FUR (Ferric Uptake Regulator) protein- in an in vivo system, the bacteria provides the FUR. Although it means we can easily increase the FUR concentration, FUR is very expensive (£130/5μg) [11].

With this in mind, these are the aims of our model:

(1) Prove that our construct would or wouldn’t work an in in vitro environment in a cell-free system.

If the system does work, we then aim to:

(1a) Find an appropriate FUR concentration for an effective efficient test, but while minimising its cost.

(1b) Find an appropriate time to take a measurement.

(1c) Show a relationship between initial Iron (II) ion concentration and amilCP concentration on a graph.

Model formation

We used mass action kinetics to model our reactions as a function of the velocity of the reaction. We first needed to establish what our reactions were. They can be seen in figure 1.

</br>

<img src="">

(1)Dimerisation of the FUR protein.

</br>

<img src="">

(2) Transcription of LacI that is inhibited by the activated FUR.

</br>

<img src="">

(3) Translation of LacI protein.

</br>

<img src="">

(4) Transcription of amiCP that is inhibited by the LacI protein.

</br>

<img src="">

(5) Translation of amilCP protein.

</br>

If we let initial FUR concentration be a constant, then more iron leads to more activated FUR, less LacI, so more amilCP. In other words, because we have two inhibitors, a lower iron concentration in saliva leads to a less intense blue, while a high iron concentration leads to a more intense blue.

We can break down our model into 9 different reactions to consider separately, then put them together once they are all defined.

For the transcription of LacI and amilCP, the transcription factors (FURactivated and LacI respectively) are not constant. We therefore cannot describe their rate’s of formation as a constant- we must consider the concentration of the transcription factors (TF). To do this, we use the hill equation, Where θ is the fraction of transcription factors that are bound to the promoter, [TF] is the concentration of the transcription factor, and Kd is the dissociation constant of the transcription factor to the promoter.

   <img src="">
   <img src="">

h is the hill coefficient (we will define it has “h” from now onwards). This can be approximated to the number of TF that binds to one promotor.

Chandran et al. discussed in their review “Mathematical modelling and synthetic biology” [1] that the rate of transcription, also known as the velocity (v) of a reaction, will be proportional to this fraction(1) if the TF is an activator, v = kTX(θ), (2) if the TF is a repressor, v = kTX(1-θ), where KTX is the maximum transcription rate of the mRNA.

   <img src="">
   <img src="">

Because both of our transcription reactions are regulated by a repressor TF, both our velocities will be equal to kTX(1-θ).

The hill equation is actually used to describe the fraction of a macromolecule saturated by ligand as a function of the ligand concentration. In this way, we can also describe our FURactivated reaction by using the hill equation, where [TF] becomes the iron (II) ions and KTX becomes initial FUR concentrations (FURt)

We can now consider each reaction. We will use v as an arbitrary variable, simply there to denote a rate of reaction. The r0x notation will be used later in the programming to remove clutter from the code

   <img src="" class="left">    

(0) Dimerisation of the FUR protein

</br>

For the FURactivated reaction, we make a quasi-steady state assumption. This is used when a part of a system reacts much quicker than the other parts, and thus reaches completion on a much smaller time scale. While FURactivated reaction will reach equilibrium in a the scale of seconds and minutes, translation and transcription takes hours, so instead of considering the velocity of the reaction, we just assume it happens instantly.

(1) Transcription of LacI mRNA (r01)

   <img src="" class="left"></br>

(2) Degradation of LacI mRNA (r02)

   <img src=""></br></br>

Note that this rate will be negative. This is because it will be decreasing the concentration of the the mRNA. We can now put them together, to form an ordinary differential equation (ODE) for the rate of mRNA production.

   <img src=""></br></br>

(3) Translation of LacI protein (r03)

   <img src=""></br>

The translation reactions are fairly simple, as they are just affected by the concentration of the mRNA of the corresponding protein and can be modelled as follows:

(4) Degradation of LacI protein (r04)

   <img src=""></br>

These two together form the ODE for the rate of LacI production:

   <img src=""></br>

(5) Transcription of amilCP mRNA (r05)

   <img src=""> </br></br> 

(6) Degradation of amilCP mRNA (r06)

   <img src=""> </br>  

These come together to form an ODE for the rate of amilCP mRNA production:

   <img src="">   </br></br></br>

(7) Translation of amilCP protein (r07)

   <img src=""></br>    

(8) Degradation of amilCP protein (r08)

   <img src=""></br>    

These two form the ODE for amilCP concentration:

   <img src=""></br> </br>   

Altogether, we have 5 reaction equations:

   <img src="" class="left"></br>
   <img src="" class="left"></br></br>
   <img src="" class="left"></br></br>
   <img src="" class="left"></br></br>
   <img src="" class="left">

Defining our constants

Now that we have our equations, we must now define our constants. This step is very important, so a lot of effort was put into ensuring we had accurate values based on the literature we could find online. These are the values we found, along with their source (note that if multiple values or a range of values were found, we indicated the value we used in a bracket {x} after the values): The rate of degradation (half-life of 60 min or faster) and the total amount of protein undergoing degradation (2 to 7%) was the same during growth and during various kinds of starvation.

For the transcription and translation rates, we used the table from Eyal Karzbrun et al [7], as it describes the rate in a cell-free system using RNA Polymerase. Transcription is described in terms of rNTP/s (1.0±0.5 -we used 1.5, as it was closer to other literature values [10]). We could then work out kTX by dividing the number of rNTPs in each mRNA by 1.5s ((1128/1.5) = 752s for LacI, (669/1.5) = 446s for amilCP) to get the time take to transcribe one mRNA, then finding it’s reciprocal ((1/752)~1.33e-3/s for LacI, (1/446) = 2.24e-3/s) to find the transcription rate.

We then did the same process for translation. Eyal Karzbrun et al [7] also describes the rate of protein synthesis in terms of a rate of amino acid production. It describes the rate as >4, but from other literature [10], their maximum is 21, so we calculated a range of values and estimated that the rate would be somewhere in the middle. LacI is 360 amino acids long [13], so 0.014≤ kTL ≤ 0.058. We used 0.04 as our value. AmilCP is 223 amino acids long [14], so 0.022≤ kTL ≤ 0.094. We used 0.06 as our value.

Solving the equations

There are fundamentally two ways to solve a series of ODEs: analytically and by using software. We first attempted to solve analytically. This involves integrating the equations so that you get a variable as a function of time. We actually have 4 variables: FURt, Fe2+. time and amilCP. The first 3 are all independent and amilCP is dependent. Trying to solve this with 4 variables is near impossible, so we let FURt and Fe2+ be constant and let time and amilCP be our two variables. In this way, after integrating fully, we would get amilCP concentration as a function of time. An attempt was made to do this. However, it proved to be unsolvable. The final equation was a second order differential of amilCP as a function of time and amilCP mRNA and even Wolfram Alpha couldn’t solve it. This being the case, we then opted to try using software, namely MatLab. Note that the k_tx and k_tl values are 0 up until a specified time. This is because in Eyal Karzbrun et al [7], they observed a delay in mRNA and protein synthesis: “Under TX/ degredation, mRNA appeared after a delay of τ0 = 15 ± 5 … We attribute the delay to the time required for TX of a single mRNA” [Fig. 2(c)]: “Assuming TX and TL are concurrent [20], we attribute the additional delay τf ≤ 5 min to the protein folding time, which is similar to previous measurements [21,22]” [Fig. 2(d)]

   <img src="" width="20%" class="left">

(a) Information flow from DNA (D) to mRNA (m) to protein (p) is carried by synthesis enzymes RNA polymerase Rp and ribosome R. Degradation by RNAse Xm and protease Xp.

(b) Exponential degradation of 200 nM mRNA with lifetime m ¼ 12 min .

(c) mRNA dynamics in the presence of 30 nM DNA exhibiting a delay and exponential accumulation with a time scale m.

(d) Kinetics of protein synthesis at various DNA concentrations.

(e) Normalized protein synthesis rate p_ syn=p_ max fitted to an exponential rise with a time scale m.

(f) Protein maximal accumulation rate as a function of DNA fitted with a MM curve using Rp ¼ 29 nM and KTX ¼ 1:5 nM.

These delays are not accounted for in our equations, so instead, we let k_tx and k_tl equal 0 up until the delay has ended- we estimated these to be 15 minutes for k_tx and 17 minutes for k_tl.

Because it is a cell-free system, transcription is 10 to 100 times slower than in vivo [7], which could be an issue, as the protein may reach steady state at a much lower concentration than we need. However, it should also be noted that because it is cell-free, we do not have to account for dilution of the proteins. Since degradation of proteins is so slow [23, 24], we can essentially ignore the degradation rates of the proteins. In our code, we set d_LacI and d_CP to 0, but we could just as well have removed r04 and r08 from the ODEs.

Unfortunately, when we ran the code, we noticed that the model had very little sensitivity to changes in FURt and iron (II) ion concentrations, as can be seen below:

Fe=0, FURt=0

</br>

   <img src="" class="left" width="48%">
   <img src="" width="48%"></br></br></br></br></br></br></br></br></br></br>

Fe=1000, FURt=0

</br>

   <img src="" class="left" width="48%">
   <img src="" width="48%"></br></br></br></br></br></br></br></br></br></br>

Fe=0, FURt=1000

</br>

   <img src="" class="left" width="48%">
   <img src="" width="48%"></br></br></br></br></br></br></br></br></br></br>

Fe=1000, FURt=1000

</br>

   <img src="" class="left" width="48%">
   <img src="" width="48%"></br></br></br></br></br></br></br></br></br></br>

As you can see, there is very little difference in changes in FURt and iron (II) ion concentrations. We tested many ranges and values, but nothing gave any sensitivity unless it was unrealistically extreme (and certainly not in the range of values we could use for FURt and we would test for Fe).

However, it did point something out to us that we hadn’t considered before: the RNA polymerase doesn’t transcribe plasmid 1 then plasmid 2- it does them both at the same time. This is shown by the amilCP mRNA concentration graph. Because there is no LacI present, there is nothing inhibiting the mRNA being transcribed so it is transcribed very quickly. Then, as LacI concentration builds up, it inhibits amilCP mRNA transcription more and more until transcription is equal to degradation (turning point). As transcription continues to decrease, degradation does not, so there is a net decrease in amilCP mRNA concentration. There is then a balance between transcription decreasing as LacI increases, while the total amount of amilCP mRNA that degrades decreases, leading to a curve as shown.

Although our model is testing in vitro values, it can be said that the lack of sensitivity is shown in our in vivo lab work too, suggesting it isn’t necessarily a mistake in the model, but perhaps an inevitability with our construct. That being said, we contemplated the idea that there was something wrong with the model, so we checked back over the constants to confirm their values. After clarifying them, we decided to try a different method altogether.

We used the Euler method instead to see if that made a difference. Instead of attempting to solve the equations, MatLab instead uses an iterative method to find each value of the amilCP after a small time period dt, then repeat up until a specified t_max This is the code we made to represent this model:

*code stuff*

This time, we decided to plot both the minimum iron (II) ion concentration (2.6e+6) and the maximum on the same graph (61.8e+6) so the difference is more clear. The blue line represents the lower bound, red represents the upper bound. This is the result:

FURt=1, dt=10, t_max=144+3600 (6 days)

   <img src="" class="left" width="100%">

As you can see (or more accurately, can’t see), there is no distinguishable difference and the blue line isn’t even visible. We checked that this wasn’t just an error by pausing the code midway through the run- MatLab plots everything it’s covered so far. As you can see, the red line is effectively in the process of covering up the blue line. This is because MatLab plots whatever is chronologically first in the code. It then plots whatever is next and importantly it appears on top of any previous plots. For us, the blue line is the lower bound, red is the upper. MatLab first plots the blue, then the red, so the red appears to cover up the blue.

   <img src="" class="left" width="100%">

We then tried to vary our FURt value. This is the result:

FURt=1000, dt=10, t_max=144*3600 (6 days)

There is an ever so slight increase in final amilCP concentration, but again, there is no distinguishable difference between the two bounds for iron levels. Finally, we put the iron (II) ion concentration values to the extremes to see what would happen. These are the results:

FURt=1000nM,Fe = 0 (blue), Fe = 1e10nM (red), dt=10, t_max=144+3600 (6 days)

   <img src="" class="left" width="100%">

As you can see, there is still no difference, despite the range of values being so extreme. For clarity, we also coded to make MatLab draw the concentration of amilCP against iron concentration after a given time t (4 hours), so see if there was any difference that we couldn’t see. This is the code we used:

bleh

These are the results

FURt=0nM

   <img src="" class="left" width="100%">

As you can see, when FURt=0, iron has no effect on amilCP concentrations, which is what we would expect.

FURt =1,000nM

   <img src="" class="left" width="100%">     

There is some variation, but it is over a very small region and it is very steep. It’s almost like a threshold iron level: below this value, there will be 777nM of amilCP, but above it there will be ~778.46nM.

FURt=1e9nM

   <img src="" class="left" width="100%">    

Interestingly, more FURt actually made the test have a smaller region on increasing (steeper increase), so we opted to try less FURt instead. We tried may values, just to see if we could see a pattern. However, less FURt just meant there was a smaller increase form 777nM. It did not widen the region to make the test work for a wider range of values. Here’s the graph showing what it looked like when FURt=1nM

FURt=1nM

   <img src="" class="left" width="100%">

We wanted to zoom in on the increasing region, so instead of testing up to 61.8e6nM, we kept testing smaller and smaller regions until we isolated just the increasing section. This graph represents iron levels varying form 0 to 0.8e6 (note that our lowest saliva iron levels is 2.6e6)

FURt=1000

   <img src="" class="left" width="100%">

Unfortunately, we couldn’t find a source showing an image of amilCP at a set concentration. This would give us an indication of when it became visible (eg, is 777nM invisible or very intense). We took this into account and tested to see if there would be any steady increase in amilCP levels on smaller time scales. We used t=1050 (30 seconds after translation begins) and t=1200 (arbitrary increase). These are the results:

time = 1050s

   <img src="" class="left" width="100%">

time = 1200

   <img src="" class="left" width="100%">

As you can see, there is still the notable threshold iron level, below which amilCP is lower, above which amilCP is higher. There is no steady increase, and certainly not over the range of values we would need to test. We simulated for many different times, and even varied FURt to see if that made a difference, but the pattern continued.

Finally, we simulated a larger time scale of 24 hours and varied FURt again, but again, the same pattern emerged:

   <img src="" width="100%">

Analysis

Summary

We set out to prove whether our construct would work in an in vitro environment in a cell-free system. We first successfully modelled our reactions in terms of ODEs. We then attempt to solve the equations in with three different methods: analytically by integration, using the MatLab ODE function and using the Euler method. Unfortunately, it was not possible to solve analytically. but both the ODE function and Euler method worked and gave similar results. The most conclusive simulations being the iron vs amilCP graphs, as they show exactly how big the difference is in amilCP and iron varies, and with different FURt and time values.

Conclusion

Our cell-free system is only sensitive over a very small region of iron levels (0 to 8e6nM) which is below our minimum iron level of 2.6e6nM. This is always the case for all FURt values and time values we tested. The model therefore quite conclusively states that the system will not work. Although our wet-lab work was in vivo, it should be noted that a similar issue was observed- that the test had no sensitivity to iron levels. We would therefore happily say that there is something that makes the system not viable. However, we we could do wet-lab work to confirm our findings.

We would create our system on a cell-free paper and vary FURt, iron and time t to the extremes, then compare with our model. It would predict that no matter the value of these three variables, the test is still not sensitive over the iron region we need it to be.

We have proved that the system does not work, and so we cannot accomplish aims 1a-1c. However, our model has perhaps highlighted another issue. It may be the case the cell-free systems are simply bad diagnostics tools. This does not absolutely conclude it, as our wet-lab work suggests it is not the cell-free nature that makes it fail, but the system itself. But there are certainly differences between the cell-free model and the wet-lab work. It may therefore be the case that even if the wet-lab predicted a successful cell-free test, it still would not have worked, as cell-free systems may just be bad a diagnostics tool.

Evaluation

As with any good model, we must evaluate the limitations of our model, highlighting any factors we haven’t accounted for and any errors that may have occurred, as well as how they could effect our conclusion.

We assumed protein degradation rate was zero: Protein degradation is very slow [23, 24], so we felt it was valid to assume it was zero, but nonetheless we could include it.

Did not account for reaction conditions: The reaction conditions of the system were ignored- temperature, pressure and pH

Other limitations on growth were not taken into account: Although we considered degradation, there are also other factors that restrict the velocity of our reaction. One example would be space. We essentially assumed there was an excess of space. However, this may not be the case, especially as we are dealing with concentrations on a μM scale. Another would be the concentrations of the enzymes were not considered. However, because we did not account for the limitations on growth, our values for amilCP are most likely higher than they would have been and even then they would be barely visible. That being the case, our model gave the reactions the best possible conditions to succeed, but the concentrations still fell short, so these assumptions are not so important- transcription is just too slow in cell-free systems.

Did not account for diffusion rates: We assumed that everything was perfectly and evenly spread among the reactants and ignored the idea that it may take time for the various parts to reach each other before they react. However, because we were considering very long periods of time and a very small surface area, this is a valid assumption and would have little effect on the outcome of our model.

Possible overuse of sources: We used Eyal Karzbrun et al [7] quite a lot in our model. This was primarily because it was a study of protein synthesis in a cell-free system, which is exactly what our project is. We also considered other sources for each of our values and took these into account when we decided on our final values for our constants.

Did not check extensively the reliability of sources: Little time was spent evaluating the reliability of sources. Although we checked that the values we found were the values we were trying to find, we did not necessarily check whether the source was particularly reliable.

Only 1 member of the team was working on the model: By only having one person working on the project, we were unable to have a second opinion at every step. In this way, mistakes are significantly more likely to occur and having multiple team members working on the model would make it much more reliable.

We started quite late, meaning there was a lot of time pressure: Because our team had no experience in modelling, it wasn’t until the UK team meet-up we attended that we got started on the model. This was when Chun Ngai Au kickstarted the modelling for us. All the work was therefore done in the course of about one month, by a highly inexperienced student. Although this couldn’t have been avoided, it should still be stated.

With these limitations on our model, it would be improper to suggest our model is 100% accurate and that our conclusion is correct. However, we feel as though we have tried to account for as much as possible, and the parts we did not account for would not have a big effect on the model. We therefore stand by our conclusion that our construct would not work in a cell-free system. However, the idea that cell-free systems are a bad way to conduct diagnostic tests is a much more significant conclusion to put forward, and especially as we only made a model for a specific case. However, we will say this may be the case.

We wish to thank George Smith and Adam Jones from the City of London School for Boys iGEM team and Chun Ngai Au from the Oxford iGEM team for guidance, assistance and discussion for our model.

References

[1] Chandran, D. et al. “Mathematical Modeling and Synthetic Biology.” Drug discovery today. Disease models 5.4 (2008): 299–309

[2] Hamed, Mazen Y., J. B. Neilands, and V. Huynh. "Binding of the ferric uptake regulation repressor protein (Fur) to Mn (II), Fe (II), Co (II), and Cu (II) ions as co-repressors: Electronic absorption, equilibrium, and 57 Fe Mössbauer studies." Journal of inorganic biochemistry 50.3 (1993): 193-210

[3] Anne Bagg and J. B. Neilands “Molecular Mechanism of Regulation of Siderophore-Mediated Iron Assimilation.” Microbiological Reviews. Dec. 1987, p. 509-518.

[4] Becskei A and Serrano L "Genetic network driven control of PHBV copolymer composition", J Biotechnol 122(1):99-121, 2006

[5] Deng, Zengqin, et al. "Mechanistic insights into metal ion activation and operator recognition by the ferric uptake regulator." Nature communications 6 (2015).

[6] Basu S et al "A synthetic multicellular system for programmed pattern formation", Nature 434:1130-1134, 2005

[7] Eyal Karzbrun et al “Coarse-Grained Dynamics of Protein Synthesis in a Cell-Free System” Physical Review Letters, Jan 28 2011, PACS numbers: 87.18.Vf

[8] Semsey et al “The effect of LacI autoregulation on the performance of the lactose utilization system in Escherichia coli” 2013 Jul; 41(13): 6381–6390

[9]Selinger, D. W. "Global RNA Half-Life Analysis In Escherichia Coli Reveals Positional Patterns Of Transcript Degradation". Genome Research 13.2 (2003): 216-223. Web. 18 Oct. 2016.

[10] H. Bremer and P. P. Dennis, in Escherichia Coli and Salmonella: Cellular and Molecular Biology, edited by F. C. Neidhardt et al. (ASM Press, Washington, DC, 1996), pp. 1553–1569, 2nd ed.

[11]https://www.abbexa.com/ecoli-ferric-uptake-regulator-protein-recombinantgclid=Cj0KCQjw95vPBRDVARIsAKvPd3J1bjFD-mag6B0xR1MofivOa6KgOigQoMOsmw8TbRZZDVBj-_C8WjAaApo5EALw_wcB

[12] Setty Y et al "Detailed map of a cis-regulatory input function", P Natl Acad Sci USA 100(13):7702-7707, 2003

[13] http://www.nature.com/nature/journal/v274/n5673/abs/274765a0.html

[14] http://parts.igem.org/Part:BBa_K592009

[15] http://www.uniprot.org/uniprot/P03023

[16] http://parts.igem.org/Part:BBa_K1163003

[18] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3711431/#gkt351-B31

[19] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3435124/table/t1-mjhid-4-1-e2012051/

[20] O.L. Miller, Jr., B. A. Hamkalo, and C. A. Thomas, Jr. “Visualization of Bacterial Genes in Action” Science 169, 392, (1970)

[21] J. Shin abd V. Noireaux, J. Biol. “Study of messenger RNA inactivation and protein degradation in an Escherichia coli cell-free expression system” Eng. 4, 9 (2010)

[22] J. A. Megerle, G. Fritz, U. Gerland, K. Jung, and J. O. Ra¨dler, Biophys. J. 95, 2103 (2008)

[23] https://slack-redir.net/link?url=http%3A%2F%2Fpubs.acs.org%2Fdoi%2Fabs%2F10.1021%2Fja954077c

[24] http://pubs.acs.org/doi/abs/10.1021/ja00230a041