Team:Athens/Model
FOR THE FULL EXPERIENCE PLEASE USE A PC :)
Model
Goal of the model
We built a kinetic model capable of monitoring dNTPs' concentrations of an E. Coli cell and the effects of overexpression of RNR and TSase enzymes. Due to scarce experimental data from the literature, the model uses statistical parameter sampling to allocate values for most of the required kinetic constants. This way a probability distribution of the dNTPs' levels is generated giving us estimations with a reasonable confidence level and dampening some of the uncertainty that the model contains. The main observations of the model are that the proposed system increases dNTPs' levels fast and in an efficient manner, RNR and TSase can be regulated by the same promoter and lastly, that the initial concentrations of dNTPs play a crucial role for the final quantities observed.
Background
RNR is the enzyme responsible for the catalysis of the conversion of the four ribonucleotides triphosphates (NTPs) into their corresponding dNTPs [1]. The mechanism of this reaction is highly complicated due to the fact that RNR is activated or inhibited through allosteric regulation by ATP, dATP, dTTP and dGTP [1]. The second enzyme, TSase plays an essential role in catalyzing the reductive methylation of deoxyuridylate (dUMP) to thymidylate (dTMP), which provides the sole intracellular de novo source of dTMP [3].
Through extensive literature research we could only find one documented model on the kinetics of RNR and TSase by Jackson et al. [2], in which the above mentioned complicated activity is described in the form of differential equations. The proposed differential equations are presented below:
$$ \begin{array}{ccl} \frac {\mathrm{d}\left( {{\mathrm{[dATP]}}} \right) } {\mathrm{d}{t} } \; &=& \; {\left( { {\frac{\frac { {{\mathrm{kcat\_ADP}} \, \cdot \, {\mathrm{[RNR]}} } \, \cdot \, {\mathrm{[ADP]}} } {{\mathrm{Km}}_{\mathrm{(ADP\_RNR)}} } } {{1} \, + \, \frac{\mathrm{[ADP]}}{{\mathrm{Km}}_{\mathrm{(ADP\_RNR)}} } } \, \cdot \, \frac {\frac{\mathrm{[dGTP]}}{\mathrm{Ka\_dGTP}} \, + \, \frac {{0.35} \, \cdot \, {\mathrm{[dTTP]}} } {\mathrm{Ka\_dTTP}} } { {{1} \, + \, \frac{\mathrm{[dGTP]}}{\mathrm{Ka\_dGTP}} } \, + \, \frac{\mathrm{[dTTP]}}{\mathrm{Ka\_dTTP}} } } \, \cdot \, \frac{1} {{1} \, + \, \frac{\mathrm{[dATP]}}{\mathrm{Ki\_dATP}} } } \right) } \; && \; {{(1)}}\\ && \\ \frac {\mathrm{d}\left( {{\mathrm{[dCTP]}}} \right) } {\mathrm{d}{t} } \; &=& \; {\left( { { { {\frac{\frac { {{\mathrm{kcat\_CDP}} \, \cdot \, {\mathrm{[RNR]}} } \, \cdot \, {\mathrm{[CDP]}} } {{\mathrm{Km}}_{\mathrm{(CDP\_RNR)}} } } {{1} \, + \, \frac{\mathrm{[CDP]}}{{\mathrm{Km}}_{\mathrm{(CDP\_RNR)}} } } \, \cdot \, \frac{\frac{\mathrm{[ATP]}}{\mathrm{ka\_ATP}} } {{1} \, + \, \frac{\mathrm{[ATP]}}{\mathrm{ka\_ATP}} } } \, \cdot \, \frac{1} {{1} \, + \, \frac{\mathrm{[dATP]}}{\mathrm{Ki\_C\_A}} } } \, \cdot \, \frac{1} {{1} \, + \, \frac{\mathrm{[dTTP]}}{\mathrm{Ki\_C\_T}} } } \, \cdot \, \frac{1} {{1} \, + \, \frac{\mathrm{[dGTP]}}{\mathrm{Ki\_C\_G}} } } \right) } \; && \; {{(2)}}\\ && \\ \frac {\mathrm{d}\left( {{\mathrm{[dGTP]}}} \right) } {\mathrm{d}{t} } \; &=& \; {\left( { { {\frac{\frac { {{\mathrm{kcat\_GDP}} \, \cdot \, {\mathrm{[RNR]}} } \, \cdot \, {\mathrm{[GDP]}} } {{\mathrm{Km}}_{\mathrm{(GDP\_RNR)}} } } {{1} \, + \, \frac{\mathrm{[GDP]}}{{\mathrm{Km}}_{\mathrm{(GDP\_RNR)}} } } \, \cdot \, \frac{\frac{\mathrm{[dTTP]}}{\mathrm{Ka\_dTTP}} } { {{1} \, + \, \frac{\mathrm{[dGTP]}}{\mathrm{Ki\_G\_G}} } \, + \, \frac{\mathrm{[dTTP]}}{\mathrm{Ka\_dTTP}} } } \, \cdot \, \frac{1} {{1} \, + \, \frac{\mathrm{[dATP]}}{\mathrm{Ki\_G\_A}} } } \, \cdot \, \frac{1} {{1} \, + \, \frac{\mathrm{[dTTP]}}{\mathrm{Ki\_G\_T}} } } \right) } \; && \; {{(3)}}\\ && \\ \frac {\mathrm{d}\left( {{\mathrm{[dTTP]}}} \right) } {\mathrm{d}{t} } \; &=& \; {\left(\frac{\frac { {{\mathrm{kcat\_TSase}} \, \cdot \, {\mathrm{[TSase]}} } \, \cdot \, {\mathrm{[dUMP]}} } {{\mathrm{Km}}_{\mathrm{(dUMP\_TSase)}} } } {{1} \, + \, \frac{\mathrm{[dUMP]}}{{\mathrm{Km}}_{\mathrm{(dUMP\_TSase)}} } } \right) } \; && \; {{(4)}}\\ && \\ \frac {\mathrm{d}\left( {{\mathrm{[dUMP]}}} \right) } {\mathrm{d}{t} } \; &=& \; {\left( { { { {\frac{\frac { {{\mathrm{kcat\_UDP}} \, \cdot \, {\mathrm{[RNR]}} } \, \cdot \, {\mathrm{[UDP]}} } {{\mathrm{Km}}_{\mathrm{(UDP\_RNR)}} } } {{1} \, + \, \frac{\mathrm{[UDP]}}{{\mathrm{Km}}_{\mathrm{(UDP\_RNR)}} } } \, \cdot \, \frac{\frac{\mathrm{[ATP]}}{\mathrm{ka\_ATP}} } {{1} \, + \, \frac{\mathrm{[ATP]}}{\mathrm{ka\_ATP}} } } \, \cdot \, \frac{1} {{1} \, + \, \frac{\mathrm{[dATP]}}{\mathrm{Ki\_C\_A}} } } \, \cdot \, \frac{1} {{1} \, + \, \frac{\mathrm{[dTTP]}}{\mathrm{Ki\_C\_T}} } } \, \cdot \, \frac{1} {{1} \, + \, \frac{\mathrm{[dGTP]}}{\mathrm{Ki\_C\_G}} } } \right) } \\ && \\ \; && \; { \, - \, \left(\frac{\frac { {{\mathrm{kcat\_TSase}} \, \cdot \, {\mathrm{[TSase]}} } \, \cdot \, {\mathrm{[dUMP]}} } {{\mathrm{Km}}_{\mathrm{(dUMP\_TSase)}} } } {{1} \, + \, \frac{\mathrm{[dUMP]}}{{\mathrm{Km}}_{\mathrm{(dUMP\_TSase)}} } } \right) } \; && \; {{(5)}}\\ && \\ \end{array} $$
The symbols and their description are listed in the Appendix. Equations (1), (2), (3), (5) contain the kinetics of the activity of RNR. Their general form contains a Michaelis-Menten term along with terms that quantify the inhibition or activation of the enzyme by the different products. To qualitatively confirm that the differential equations from Jackson et al correctly describe the allosteric regulation of RNR, we referred to a paper by Henderson et al, in which the influence of different effectors on the activity and specificity of RNR is documented [4]. Lastly equation (4) represents the activity of TSase in a regular Michaelis-Menten form. It is worth noting that the original equation considers the kinetic parameters of equation (4) as dependent on the concentrations of other biochemical species, but at the present model are assumed constant due to the marginal quantities of those species in constant with the overproduced dNTPs.
In addition, in order to describe the expression of RNR and TSase we included two more equations. These equations consider the kinetics of association-disassociation of the inducer and the transcription of the TU to mRNA, in a steady state due to the difference in time scales (these reactions occur in the scale of seconds while the system operates in the scale of hours). The equations that were constructed are presented below:
$$ \frac{d R N R}{d t}=a_{R N R} \frac{I P T G}{K_{d, L c c}+I P T G}-d_{2, R N R} R N R \quad {{(6)}} $$ $$ \frac{d \text { TSase }}{d t}=a_{ \text { TSase }} \frac{\text { AraC }}{K_{\text {d, trac }}+\text { AraC }}-d_{2,\text { TSase }} \text { TSase } \quad {{(7)}} $$ $$ a_{R N R}=k_{2, R N R} \frac{k_{1, rnr}}{d_{1, rnr}} C_{N} \quad K_{d, \mathrm{Lac}}=k_{o f f, \mathrm{Lac}} / k_{o n, \mathrm{Lac}} \quad a_{T S a s e}=k_{2, T \text { Sase }} \frac{k_{1, t h y A}}{d_{1, t h y A}} C_{N} \quad K_{d, A r a C}=k_{o f f} / k_{o n} $$
Equation (6) was based on documentation of the T7-LacO promoter from the iGEM team Edinburgh UG [5] with some minor modifications. In this equation, k1 is proportional to the copy number of the plasmid, the concentrations of the genes is mostly considered constant and depends on the origin of replication and the copy number of the plasmid, RNA polymerase and ribosomes are many enough to not hinder the kinetics and the binding and unbinding is assumed not to play a role in the kinetics. Equation (7) was based on the paper by Schleif et al [6] and uses an AraC inducer instead of IPTG.
All the above equations are based on the intracellular dynamics of the system. In order to have estimations for a liquid bacterial culture, we used Monod equations to upscale the results. An additional negative term -pX is used, whose value would be determined by fitting to experimental data provided by wet lab experiments. This term aims to simplistically quantify the anticipated effects of mutagenesis that the activity of elevated concentrations of RNR and TSase will have on the bacterial concentration of the culture. These effects can be translated as loss of the plasmid, loss of function of the designed genetic circuit or even unbalanced operation of the bacteria itself. Unfortunately, we did not have wet lab data to actively use these equations so the results presented below are representative of the intracellular kinetics of the system.
$$ \begin{array}{c} \mu=\mu_{\max } \frac{S}{K_{S}+S} \\ \frac{d X}{d t}=\mu X-k_{d} X-p X \\ \frac{d S}{d t}=-\frac{1}{Y_{X / S}} \mu X \end{array} $$
Technical details of model
The above described model operates under some assumptions and simplifications that should be pointed out. Firstly, the bacterial cell is treated as an open system in a steady state meaning that some biochemical species are provided at a constant concentration (i.e. ADP, CDP, GDP, UDP and ATP) and that the enzymatic activity of RNR and TSase are considered the rate limiting step of the whole metabolic pathway. An aftereffect of this assumption is that the deoxynucleoside di- and triphosphates are readily interconvertible due to the kinases’ rapid activities. This means that for example, RNR may produce dADP from ADP but the final presented species is dATP because the phosphorylation of the species happens rapidly.
A simplification that is part of the model is that we do not take into consideration the concentrations of dNTPs coming from salvage pathways. That is because the salvage pathway is mostly used when the cell is in need and in our modeled system, the de novo production of dNTPs will oversaturate the cell.
One of the main issues of the proposed model is that it has a vast number of kinetic parameters, which do not have documented values derived from experimental results. Specifically, the reference for the RNR kinetics [2] uses a lot of unpublished data or educated guesses to derive values for the kinetic parameters and the whole model is built based on mammalian cells. As the paper states in many instances, the explored effects are more qualitative validations rather than exact descriptions of the system. Based on some experimental data, we know that the kinetic parameters for a bacterial cell differ a lot, sometimes 2-3 orders of magnitude, from the respective values used for the mammalian cell [7].
However, the literature lacks key experimental data that would replace the values of the reference. That means that the model is characterized by a lot of uncertainty, which translates to vague results and estimations. To combat this, a parameter sampling method must be used in order to create results with a variance and some confidence level. This method requires three additional steps:
- We used the protein database BRENDA [7] in order to define the parameter space. The goal was to couple each kinetic parameter with a mean value and a standard deviation (sd) value. Moreover, the values should be based on E. Coli cells. In the very likely event that no sd value was documented, then we produced one, based on the mean values documented for the same modifier (i.e. substrate, activator or inhibitor species) but for different organisms. In the less likely case that the database did not have a mean value for a parameter, then we produced one based on the mean values documented for the same modifier (i.e. substrate, activator or inhibitor species) but for different organisms.
- Next we created a parameter scan task in COPASI (with which we have built the whole of the kinetic model) that would do parameter sampling for each kinetic parameter, from their own normal distributions. The normal distributions derived from the mean and the sd value of each parameter.
- Through this framework we created a set of 100000 (100K) models for each case study we had to explore.The model ran for a simulation time of 30000s (~8 hours). It is worth noting that the number of generated models could be higher but the timeframe of the iGEM project and the computational power available (standard laptop running Windows 10 with intel i5 processor) would not allow for a lengthier simulation time.
The proposed model framework can be found in our github page. Information on the mean and sd values for the different kinetic parameters along with references to their documentation can be found in the Appendix.
Based on the presented background, assumptions and preparation of the model, the different cases we wanted to test could begin.
Figure 1. Left: dNTP concentrations with respect to time. The filler lines represent the standard deviation of each time point (dTTP showcases marginal standard deviation). Right: frequency histogram of the calculated concentrations of dATP at 30000s.
| Species |
|---|
| dATP 0.43 ± 0.10 mM |
| dCTP 0.24 ± 0.02 mM |
| dGTP 0.74 ± 0.10 mM |
| dTTP 0.25 ± 0.01 mM |
Case 1: Fixed RNR and TSase levels
The initial test that was performed was the observation of the system in fixed RNR and TSase concentrations, a state indicative of the WT E. Coli cell. As such, equations (7) and (8) were not included and 100K models were produced. The results are presented in Figure 1. For simplicity reasons, we will present analytical histograms only for dATP but will comment on the values of other dNTPs too.
dATP, dGTP concentrations reach values barely above 0.4 mM while dCTP and dTTP values are lower but at the same order of magnitude. The system is quite slow in the formation of the desired quantities of dNTPs and showcases an irregularity; dTTP levels are unexpectedly low. This can be attributed to the slow formation of dUMP which instantly turns into dTTP and subsequently showcases low values of dUMP. This irregularity persists between all the different cases performed and will be discussed in the conclusions section of the page. Another observation is that the sd values are very small meaning that despite substantial changes to the values of the kinetic parameters the system reaches similar quantities in every model run. Lastly, dGTP concentrations seem to slightly overpass dATP concentrations. This is not expected but indeed, if the model runs for longer dATP ends up overpassing dGTP’s concentrations.
Figure 2. Left: dNTP concentrations with respect to time. The filler lines represent the standard deviation of each time point (dTTP showcases marginal standard deviation). Right: frequency histogram of the calculated concentrations of dATP at 30000s.
| Species |
|---|
| dATP 219.3 ± 57.1 mM |
| dCTP 1.84 ± 0.32 mM |
| dGTP 27.2 ± 7.3 mM |
| dTTP 0.51 ± 0.06 mM |
Case 2: Transient RNR and TSase levels, with different regulatory promoters
The next step that was examined contained the equations (7) and (8) making the concentrations of TSase and RNR transient. This version of the model aims to estimate the effects of the genetic circuits that will be introduced to the E. Coli cell. The 100K models were produced and the results are presented in Figure 2.
The results show that the system greatly overproduces dNTPs in contrast with the WT E. Coli cell. Specifically, dATP concentration is very high, dCTP and dGTP is one order of magnitude lower while dTTP had the smallest change but still overpasses the desired concentration 0.4 mM while doubling in size. The sd values remain relatively low so the system can still be considered stable. The variance observed also implies that despite substantial differences in the values of the kinetic parameters the system results in a desirably confined space.
Figure 3. Left: dNTP concentrations with respect to time. The filler lines represent the standard deviation of each time point (dTTP showcases marginal standarad deviation). Right: frequency histogram of the calculated concentrations of dATP at 30000s.
| Species |
|---|
| dATP 219.3 ± 56.9 mM |
| dCTP 1.84 ± 0.32 mM |
| dGTP 27.2 ± 7.3 mM |
| dTTP 0.46 ± 0.08 mM |
Case 3: Transient RNR and TSase levels, with the same regulatory promoter
Due to changes in the wet lab design, during the later cycles of the project, a new system was built which contained one common promoter for RNR and TSase. This initiated a second case study, with which a last version of the model was built to observe how the system would react when the same promoter was used. Specifically, the T7-Lac promoter was used (IPTG inducer) and the parameters in equation (8) were changed accordingly. The results of this version of the model are shown in Figure 3.
As shown above, the system has very similar results with the previous case study. dATP levels remain high while dTTP mean value showcases a small drop, that however is between the original standard deviation. It is certain that these results are positive about the new wet lab design. However, it should be noted that the model contains some assumptions without which the results would be different. Allosteric modifiers of TSase, alternative metabolic pathways as well as limitations in the concentration of ribonucleotides could be some of the examples that the system will behave differently in the wet lab experiments. Only thing left to do is begin the experiments!
Conclusions
Based on the above observations a general set of conclusion can be derived:
- The modeled system had great changes on dNTP concentrations and rate of formations when kinetics for overexpression of RNR and TSase were introduced. The time needed for the production of the desired quantities of dNTPs was also vastly reduced, meaning that the project’s goals of efficiency and cost effective production were plausible.
- It is possible to regulate the overexpression of RNR and TSase with the same promoter, without having negative changes to the final results.
- dTTP levels were unexpectedly low across all case studies. This also resulted in low levels of dUMP which was rapidly turned into dTTP. While not much can be said about this characteristic of the system, it might be possible that the cell will experience shortages of UDP ribonucleotides, affecting its normal biological function at the same time.
- The system was very sensitive to the values of the initial concentrations of dNTPs, giving vastly different results. This should be taken into consideration when performing the wet lab measurements.
These observations come with a considerable level of uncertainty due to the scarcity of documented literature values as well as experimental data. Through measurements from our wet lab experiments, we aim to reduce this uncertainty by determining the values of some of the kinetic parameters. Here is how we would do it.
Future planning
RNR and TSase enzymes were originally linked with GFP and RFP enzymes respectively. With the help of fluorescent measurements we would better determine the kinetic parameters of the equations (6) and (7) for the overexpression of the enzymes. Additionally, dNTPs concentration measurements with the use of HPLC readings would provide levels of each dNTP at different time points of the culture and help us better calculate the kinetic parameters of equations (1) through (5).
Appendix
| Kinetic Parameter | Explanation |
|---|---|
| Km_ADP | Michaelis-Menten constant of RNR kinetics. The substrate is ADP. |
| Km_CDP | Michaelis-Menten constant of RNR kinetics. The substrate is CDP. |
| Km_GDP | Michaelis-Menten constant of RNR kinetics. The substrate is GDP. |
| Km_UDP | Michaelis-Menten constant of RNR kinetics. The substrate is UDP. |
| kcat_ADP | Catalytic rate constant of RNR kinetics. The substrate is ADP. |
| kcat_CDP | Catalytic rate constant of RNR kinetics. The substrate is CDP. |
| kcat_GDP | Catalytic rate constant of RNR kinetics. The substrate is GDP. |
| kcat_UDP | Catalytic rate constant of RNR kinetics. The substrate is UDP. |
| Ki_A_A | Inhibition constant of RNR kinetics. The substrate is ADP and the inhibitor is dATP. |
| Ki_G_A | Inhibition constant of RNR kinetics. The substrate is GDP and the inhibitor is dATP. |
| Ki_G_T | Inhibition constant of RNR kinetics. The substrate is GDP and the inhibitor is dTTP. |
| Ki_G_G | Inhibition constant of RNR kinetics. The substrate is GDP and the inhibitor is dGTP. |
| Ki_C_A | Inhibition constant of RNR kinetics. The substrate is CDP and the inhibitor is dATP. |
| Ki_C_T | Inhibition constant of RNR kinetics. The substrate is CDP and the inhibitor is dTTP. |
| Ki_C_G | Inhibition constant of RNR kinetics. The substrate is CDP and the inhibitor is dGTP. |
| Ka_dGTP | Activation constant of RNR kinetics. The activator is dGTP. |
| Ka_dTTP | Activation constant of RNR kinetics. The activator is dTTP. |
| Ka_ATP | Activation constant of RNR kinetics. The activator is ATP. |
| kcat_TSase | Catalytic rate constant of TSase. The substrate is dUMP. |
| Km_TSase | Michaelis-Menten constant of TSase. The substrate is dUMP. |
| Kinetic Parameter | Mean Value | Standard Deviation | Units | Reference |
|---|---|---|---|---|
| Km_ADP | 0,07 | 0,011 | mM | * |
| Km_CDP | 0,05 | 0,008 | mM | * |
| Km_GDP | 0,16 | 0,017 | mM | * |
| Km_UDP | 0,22 | 0,054 | mM | * |
| kcat_ADP | 1,4 | 0,3 | 1/s | [8] [2] |
| kcat_CDP | 8,0 | 0,3 | 1/s | [8] [2] |
| kcat_GDP | 3,0 | 0,3 | 1/s | [8] [2] |
| kcat_UDP | 4,4 | 0,3 | 1/s | [8] [2] |
| Ki_A_A | 0,0004 | 0,0002 | mM | * |
| Ki_G_A | 0,004 | 0,0002 | mM | * |
| Ki_G_T | 1,4 | 0,0005 | mM | * |
| Ki_G_G | 0,0098 | 0,0006 | mM | * |
| Ki_C_A | 0,01 | 0,0002 | mM | * |
| Ki_C_T | 0,05 | 0,0005 | mM | * |
| Ki_C_G | 0,006 | 0,0006 | mM | * |
| Ka_dGTP | 0,0005 | - | mM | [2] |
| Ka_dTTP | 0,04 | - | mM | [2] |
| Ka_ATP | 0,4 | - | mM | [2] |
| kcat_TSase | 1,32 | 0,02 | 1/s | [9] |
| Km_TSase | 0,5 | 0,2 | uM | [9] |
(*): an excel file with database values and comments for every uskinetic parameter can be found here. The code to reproduce our results can be found here.
1. Torrents, E. (2014). Ribonucleotide reductases: essential enzymes for bacterial life. Frontiers in Cellular and Infection Microbiology, 4. https://doi.org/10.3389/fcimb.2014.00052
2. Jackson, R. C. (1984). A kinetic model of regulation of the deoxyribonucleoside triphosphate pool composition. Pharmacology & Therapeutics, 24(2), 279–301. https://doi.org/10.1016/0163-7258(84)90038-x
3. Liu, J., Schmitz, J. C., Lin, X., Tai, N., Yan, W., Farrell, M., Bailly, M., Chen, T., & Chu, E. (2002). Thymidylate synthase as a translational regulator of cellular gene expression. Biochimica et Biophysica Acta (BBA) - Molecular Basis of Disease, 1587(2-3), 174–182. https://doi.org/10.1016/s0925-4439(02)00080-7
4. Henderson, J. F., & Paterson, A. R. P. (1973). ENZYMATIC REDUCTION OF RIBONUCLEOTIDES. Nucleotide Metabolism, 244–263. https://doi.org/10.1016/b978-0-12-340550-0.50024-1
5. https://2017.igem.org/Team:Edinburgh_UG/Model
6. Schleif, R. (2010). AraC protein, regulation of the l-arabinose operon inEscherichia coli, and the light switch mechanism of AraC action. FEMS Microbiology Reviews, 34(5), 779–796. https://doi.org/10.1111/j.1574-6976.2010.00226.x
7. https://www.brenda-enzymes.org/index.php
8. Ravichandran, K., Olshansky, L., Nocera, D. G., & Stubbe, J. (2020). Subunit Interaction Dynamics of Class Ia Ribonucleotide Reductases: In Search of a Robust Assay. Biochemistry, 59(14), 1442–1453. https://doi.org/10.1021/acs.biochem.0c00001
9. Islam, Z., Gurevic, I., Strutzenberg, T. S., Ghosh, A. K., Iqbal, T., & Kohen, A. (2018). Bacterial versus human thymidylate synthase: Kinetics and functionality. PLoS ONE, 13(5), e0196506. https://doi.org/10.1371/journal.pone.0196506