Team:RHIT/Model
Modeling
Introduction:
This model would predict the fluorescence level over time in single cells with respect to different concentrations of vitamin B9/B12 inputs. It was able to predict the trend of fluorescence level over 20 hours after inoculated cells were taken from their cultures and placed in M9 minimal media.
The system consists of a riboswitch on mRNA that would partially inhibit the production of GFP upon binding with the respective vitamin (B9/B12) (Beisel and Smolke, 2009). Each riboswitch binding vitamin B9 was found to bind two molecules (Batey et al., 2011) while the one binding B12 was found to bind only one (Batey et al., 2017). The mRNA produced by transcription would form hairpins that would stabilize upon binding with the vitamins. These hairpins would then reduce the amount of the GFP produced downstream by a factor depending on individual riboswitches. The mRNA and GFP molecules would then degrade over time. They would also diffuse due to cell growth according to (Beisel and Smolke, 2009). The production, degradation, and diffusion would result in several equilibriums. By specifying the initial amount of vitamin input, the amount of GFP level over time can be modeled.
It was simplified with several assumptions. The cell growth was simplified to a number of segments and the rate of vitamin transporting through the membrane was assumed to be instant. Since calibration curve relating OD630 measurements and cell counts was not constructed, it was difficult to estimate the concentration of vitamin in each cell experimentally. The effect of the decrease in concentration of vitamin outside of the cell due to the increase of the cell volumes was assumed zero. The pattern of production of mRNA and GFP and the growth rate of cell were also generalized to be a few constants into a few segments.
Definitions of Variables & Constants:
Functions for the Riboswitches:
R1(t) = [unfolded RNA], mM
R2(t) = [folded RNA], mM
R3(t) = [folded RNA bond to vitamin B9], mM
B(t) = [Vitamin B9 within the cell], mM
F(t) = [GFP (BBa_E0040)], mM
Constants:
Vocabulary:
kr = Production rate of mRNA, M/s^-1
k1 = Rate that mRNA folds, s^-1
k1’ = Rate that folded mRNA unfolds, s^-1
k2_B9 = Rate that folded mRNA binds with vitamin B9, M^-2 * s^-1
k2_B12 = Rate that folded mRNA binds with vitamin B12, M^-1 * s^-1
k2’ = Rate that mRNA that bond to vitamin unbinds, s-1
kA = Production rate of GFP by unfolded mRNA, s^-1
kB = Production rate of GFP by folded mRNA, s^-1
kdm = Degradation rate of mRNA, s^-1
kdP = Degradation rate of GFP, s^-1
μ = Diffusion rate inside the cell due to cell multiplying, s^-1
Constant values:
kr = 10^-9 M*s^-1
k1 = 0.1 s^-1
k1’ = 1 s^-1
k2_B9 = 34375 M^-2 * s^-1
k2_B12 = 34375 M^-1 * s^-1
k2’ = 1 s^-1
kA = 0.035 s^-1
kB = 0.3*kA s^-1
kdm = 2.31*10-3 s^-1
kdP = 5.8346*10-6 s^-1
μ = 6.7956*10-4 s^-1 (Other phases of cell growth)
μl1 = 0 s^-1 (lag1 phase of cell growth)
μl2 = -1.6554*10-5 s^-1 (lag2 phase of cell growth)
Differential Equations for vitamin B9 (THF) Riboswitches:
dR1/dt=kr-k1*R1+k1’*R2-(kdm+μ)*R1
dR2/dt=k1*R1-k2_B9*R2*B^2+k2’*R3-k1’*R2-(kdm+μ)*R2
dR3/dt=k2_B9*R2*B^2-k2’*R3-(kdm+μ)*R3
dB/dt≈0
dF/dt=kA*R1+kB*(R2+R3)-(kdP+μ)*F
Differential Equations for vitamin B12 Riboswitches:
dR1/dt=kr-k1*R1+k1’*R2-(kdm+μ)*R1
dR2/dt=k1*R1-k2_B12*R2*B+k2’*R3-K1’*R2-(kdm+μ)*R2
dR3/dt=k2_B12*R2*B-k2’*R3-(kdm+μ)*R3
dB/dt≈0
dF/dt=kA*R1+kB*(R2+R3)-(kdP+μ)*F
Results and Discussion:
While many assumptions were made, the graphs would serve as an approximation of the trend of fluorescence level over time. The mRNA molecules, produced by transcription of the DNA, would lead to the production of GFP. The mRNA molecules and GFP would also decay at the same time while the cell growth due to division would also decrease their concentration. An equilibrium of production and decay could then be reached. While the mRNA form hairpins and bind with vitamin B9/B12, the production rate of GFP would be reduced (Beisel and Smolke, 2009), shifting the equilibrium and lowering the concentration of GFP. The different intervals served to model the different phases of cell growth. In each phase the production rate of mRNA and GFP, the growth rate, and the rate of decay of several items should all differ in various degrees due to changes in the cell and the environment, causing the fluorescence to change in different phases. This model captured parts of these differences.
Figures:
Vitamin B9 riboswitches with different concentration input of vitamin B9:
Vitamin B9 riboswitches zoomed in at the point where equilibrium is almost reached:
Vitamin B12 riboswitches:
Vitamin B9 riboswitches zoomed in at the point where equilibrium is almost reached:
Several constants for specific riboswitches can only be found experimentally but the number of experiments conducted was very limited, so only one set of graphs for each of the vitamin B9 and B12 riboswitches were made.
The model predicts that the fluorescence would fall in the second interval (after the cells were transferred into the M9 minimal media, growth phase: lag1). The fluorescence would then experience a continuous increase in the subsequent interval (growth phase: lag2). Then it would dramatically fall until reaching a new equilibrium (growth phases: exponential and stationary).
It was assumed that the promoters were not active in the lag1 phase (second interval), which meant almost no mRNA or GFP was made according to (Alon et al., 2013). When the promoter came back to function in the lag2 phase (third interval) and the cells barely grew in size (Alon et al., 2013), the fluorescence largely increased. When the cell finally began to grow in the exponential phase (fourth interval), the concentration of fluorescence fell again until reaching an equilibrium.
This was a model with many assumptions that could be improved further. Whether any of the productions was reduced to ~0 in the lag1 phase was not very clear. The conducted experiment indicated that fluorescence might peak between 2-4 hr, much earlier than the end of lag2 phase predicted by the model. The concentration used for the models for the vitamin B12 riboswitches could be too high or the values of kA and kB could be too different from the actual values. Therefore, the graphs stuck too close to each other. The rate of cell division in the original medium and the stationary phase was also unclear, so the diffusion constant for the exponential phase was used for them. It should be higher in the original medium with more nutrients and lower in the stationary phase due to the decrease in nutrients. The significance of these, however, could be reduced by only looking at the fluorescence level at the end of the exponential phase, where the diversion of fluorescence caused by different concentrations of vitamins would become very obvious. The probable lower division rate and mRNA/GFP production rate in the stationary phase should have more complex results that were not able to be modeled due to lack of information.
Methods and Coding:
Calculations of values of the constants:
Degradation rate of RNA:
The half-life of a RNA molecule in exponential phase is 4.1 minutes and 7.8 minutes in the stationary phase as in (Xie et al., 2015). A basic differential equation to find out the degradation constant was constructed:
kdm =degradation constant of mRNA
t = half-life of mRNA, min
(d[RNA])/dt=-kdm*[RNA]
Solving it gives:
[RNA]=e^(-kdm*t) kdm=(ln(2))/t
Calculating half-life: (exponential phase)
kdm= (ln(2))/(4.1*60)
The degradation constant of mRNA:
kdm=0.00282 s^(-1)
It was assumed that the RNA degradation rate does not change no matter of the folding pattern.
The degradation constants was calculated with the same method:
Other degradation constants:
The average E.coli division time in M9 minimal media is about 67 minutes: (BNID 107132, Milo et al., 2010. )
µ= (ln(2))/(67*60)=1.7242*10^-4 s^-1
According to the iGEM registry, the half-life of the GFP (BBa_E0040) is about 33 hours:
kdP= (ln(2))/(33*3600)=5.8346*10^-6 s^-1
The E.coli growth in length in the lag2 phase is approximately 15% (Alon et al.), which would result in an approximately 15% increase in volume. The linearly approximated μ for the interval (~5.4 hr) is therefore:
(1+μ)^(5.4*3600)=1/(1+0.15)
μ_l2=-1.6554*10^-5 s^-1
Rate of protein production:
For the Smu THF riboswitch with a chimeric metE platform, the folded form of RNA produces about 30% of the amount the unfolded form does (Batey et al., 2011). For the L. casei one, it was suggested to be a little above 25% (Leigh, 2015).
Therefore, it was assumed that it would be:
kB=C*kA
Where C ≈ 0.3 for the B9 Riboswitch with metE expression platform. (Batey et al., 2011). It would vary significantly by riboswitches, but insufficient data were found for most of them, so 0.3 was used.
Rate constant controlling mRNA folding:
Say:
K1 = k1’/k1
K2 = k2_B9/k2’
Then the relationship between K1 and K2: (EC_50: Concentration needed for half-maximal expression response) (Beisel and Smolke, 2009)
EC_50=(1+K1*kdMA/kdMB)/K2
It was assumed the constants kdm = kdMA = kdMB, therefore the equation is reduced to:
EC_50=(1+K1)/K2
The EC50 of the Smu THF riboswitch with the chimeric metE expression platform with (6S)-THF is about 32*10-6 M. (Batey et al., 2011)
Reorganizing the equation:
K2=(1+K1)/(32*10^-6)
K2=31250*K1+31250
Since the B9 Riboswitches bind two THF molecules at once, the unit of k2_B9 would be M^-2*s^-1.
Assumed Constants:
Values of k1, k1’, k2_B9, k2_B12, and k2’:
There was not really a way to easily find out k1, k1’, k2_B9, k2_B12, and k2’. They had to be assumed from the ratios.
From our previous calculations it was assumed: K1 = 0.1, then K2 = 34375
Since binding with vitamin would stabilize the hairpin on mRNA, the reaction constant would be it was assumed that
k1 = 0.1 s-1, k1’ = 1 s^-1, k2_B9 = 34375 M^-2*s^-1, k2’ = 1 s^-1
There was very little information on the binding strength of cyanocobalamin, the type of vitamin B12 that was used in the experiment, with the vitamin B12 (AdoCbl) riboswitches which primarily bind to adenosylcobalamin. Therefore, the same numerical value was used for both k2_B9 and k2_B12 while the units are different.
Value of kA:
The translation rate of ribosome was approximated to be 8.4 aa*s-1 while Ribosome per mRNA is about 3.46 Ribosome/100 codon. (aa=amino acid) (Siwiak and Zielenkiewicz, 2013) Length of GFP (BBa_E0040) = 720pb = 240 aa*GFP^-1
Therefore:
kA = Translation rate/number of amino acids = 8.4 aa*s-1 /240 aa*GFP^-1 = 0.035 GFP*s^-1
Value of kr:
It was very difficult to find out information that correlates the promoter strength to transcription rate. It was assumed its value to be 10^-9 M*s^-1. (Beisel and Smolke, 2009)
Coding
The model was coded using Matlab.
The code was done in four time intervals using sub-functions for reiterations of differential equations.
The initial conditions of the subsequent intervals were set to be the same as the last value of the previous intervals.
In the first interval, the initial concentrations of all items (R1, R2, R3, Bin, Bout, F) were set to zero and let the system get to an equilibrium. Then, the next three time intervals were calculated sequentially. The second and third intervals represented the lag1 and lag2 phases and the fourth interval represented the exponential and stationary phases of E.coli growth. In the second interval (lag1 phase), it was assumed that no mRNA containing riboswitches was produced, no cell division occurred and the membrane pumps were not functional due to lack of ATP, therefore kr = 0 and µ = 0, and no vitamin B is present in the system. (Alon et al., 2013) Then, the vitamin was “added” as an initial condition of B of the third interval since the cell was assumed not producing any GFP at all (Alon et al., 2013). It was assumed that vitamin B enters the system almost instantly. In the third interval (lag2 phase), it was assumed that there were still no cell divisions, but they grew in size over a long period of time and mRNAs were being produced. Therefore, µ was small and kr > 0. Modeling the change of kr in this period was difficult, so it was assumed to begin at its maximum value. The length of the two intervals is about 3 and 5.4 hours in M9 minimal media. (Alon et al., 2013) In the fourth interval, it was assumed that the cell behaved the same in both exponential and stationary phases due to lack of information on many components in the stationary phase. It was assumed that the division rate of the cells did not change but the death rate caught up and the GFP production rate per cell was not affected by the dead cells.
The lengths of the intervals were 50000s (~14hr, a time sufficient for the GFP concentration in cells to reach equilibrium), 10800s (3hr), 19440 (5.4hr), and 41760 (11.6hr, just for that the time add up to 20hr, which was the experimental time).
After calculating the concentrations of GFP, the concentrations of GFP was multiplied by 79429, the amount of fluorescence (Au) per mM according to the iGEM registry.
Citations
Beisel, C.L. & Smolke, C.D. Design Principles for Riboswitch Function. PLoS. 2009.
https://doi.org/10.1371/journal.pcbi.1000363
Chen, H., Shiroguchi, K., Ge, H. & Xie, X.S. Genome-wide study of mRNA degradation and transcript elongation in Escherichia coli. Mol Syst Biol. 2015; 11: 781.
https://doi.org/10.15252/msb.20145794
Leigh, J. Rational Design, Re-engineering and Characterisation of Tetrahydrofolate Riboswitches in Bacteria (Doctoral Dissertation). University of Manchester eScholar. 2015.
https://www.escholar.manchester.ac.uk/api/datastream?publicationPid=uk-ac-man-scw:281096&datastreamId=FULL-TEXT.PDF
Madar, D., Dekel, E., Bren, A., Zimmer, A., Porat, Z., Alon, U. Promoter Activity Dynamics in the Lag Phase of Escherichia coli. BMC Syst Biol. 2013; 7: 136.
https://doi.org/10.1186/1752-0509-7-136
Milo, R., Jorgensen, P., Moran, U., Weber, G. & Springer, M. BioNumbers-the Database of Key Numbers in Molecular and Cell Biology. BNID 107132. Nucl. Acids Res. 2010; 38: D750-D753.
https://doi.org/10.1093/nar/qkp889
Polaski, J.T., Webster, S.M., Johnson, J.E. & Batey, R.T. Cobalamin riboswitches exhibit a broad range of ability to discriminate between methylcobalamin and adenosylcobalamin. JBC. 2017; 292: 11650-11658.
https://doi.org/10.1074/jbc.M117.787176
Trausch, J.J., Ceres, P., Reyes, F.E. & Batey, R.T. The structure of a tetrahydrofolate-sensing riboswitch reveals two ligand binding sites in a single aptamer. Structure. 2011; 19: 1413-1423.
https://doi.org/10.1016/j.str.2011.06.019
