DNA-Based System Models

Introduction

We created a kinetic model to simulate the dynamics of our cell free DNA-based system and to inform the design process. Using this model, we made changes to our original design, most notably the addition of an amplification step to increase anticoagulant production. We modelled the system mechanistically to make it as accurate as possible, using mass action kinetics and Michaelis-Menten kinetics to model the reactions. The goal of this model was to answer the following questions:

1. How fast is hirudin produced?
2. How leaky is the TetR promoter? What is the foldchange between a negative and positive test?
3. Difference between producing hirudin directly and producing TEV to cleave inactivated hirudin?

Methodology

To model the biochemical reactions, we formulated a system of ordinary differential equations (ODEs) that are solved in MATLAB.

The Cell-free DNA system is modelled mechanistically by using dissociation constants found in literature for TetR-TetO binding and RNAP-pTet binding. TetR dimerisation is modelled to introduce TetR monomers that may act to competitively inhibit Cruzipain in cleaving the TetR dimer. Furthermore, ribosomes are modelled to decay over time in the cell-free extract to saturate protein expression at around 2-3h. Because only dissociation constants are obtained, certain rate constants are estimated from the equilibrium constants.

The following parameters were used:

Parameter	Variable Name	Value	Reference*
Dissociation Constant of RNAP-DNA	\(K_{d,rd}\)	\(1.26*10^{-8}\)	E
Association rate of RNAP-DNA	\(k_{a,rd}\)	\(5*10^6\)	E
Dissociation rate of RNAP-DNA	\(k_{d,rd}\)	\(k_{a,rd}*K_{d,rd}=0.0632\)	E
Isomerization Rate of RNAP-DNA	\(k_{2,rd}\)	\(1.03*10^{4}\)	E
Transcription Rate of Hirudin	\(k_{TX}\)	\(0.0033\)	C
Degradation Rate of mRNA	\(k_{deg,mRNA}\)	\(log_2/(4*60)=0.0029\)
Translation Rate of Hirudin	\(k_{TL}\)	\(0.0615\)	C
Dissociatin Constant for the Ribosome and DNA	\(K_{Ribo}\)	\(10^{-9}\)	–
Degradation Rate for the Ribosome	\(k_{deg,Ribo}\)	\(7.50*10^{-5}\)	F
Association Rate for the promoter	\(k_{a,p}\)	see below	D
Dissociation Rate for the promoter	\(k_{d,p}\)	see below	D
Activation Rate for TetR	\(k_{2,TetR}\)	see below	D
Dissociation Constant for TetR-TetR and 2TetR	\(K_{d,TetRd}\)	\(10^{-8}\)	B
Association Rate for TetR-TetR and 2TetR	\(k_{a,TetRd}\)	\(10^9\)	B
Dissociation Rate for TetR-TetR and 2TetR	\(k_{d,TetRd}\)	\(K_{d,TetRd}*k_{a,TetRd}=10\)	B
Michaelis Constant for Cruzipain	\(K_{m,c}\)	\(5.8*10^{-6}\)	A
Catalysed rate of reaction for Cruzipain	\(k_{cat,c}\)	\(10.8\)	A
Association rate for Cruzipain	\(k_{a,c}\)	\(10^7\)	Estimated
Dissociation rate for Cruzipain	\(k_{d,c}\)	\(K_{m,c}*k_{a,c}-k_{cat,c}=47.20\)	A
Cleavage reaction rate for the second cleavage site in the TetR dimer	\(k_{a,c2}\)	\(k_{a,c}*10^{3}\)	A
Dissociation rate for Cruzipain	\(k_{d,c2}\)	same as \(k_{d,c}=47.20\)	A

The following species were used

Species Name	Symbol Used
DNA	\(DNA\)
mRNA	\(mRNA\)
RNA Polymerase	\(RNAP\)
A specific binding of RNA Polymerase to a promoter region to form a closed complex	\(DNA:RNAP\)
The activated RNA Polymerase binded to DNA	\(DNA:RNAPa\)
Ribosome	\(Rib\)
Hirudin	\(Hirudin\)
Tet Repressor Protein	\(TetR\)
TetR dimer	\(TetR_{2}\)
DNA binded to TetR dimer (Intermediate Complex)	\(DNA:TetR_{2}\)
DNA binded to TetR dimer (Specific Complex aka activated form)	\(DNA:TetRa_{2}\)
Cruzipain	\(Cruzipain\)
DNA binded to activated TetR dimer and Cruzipain	\(DNA:TetRa_{2}:Cruzipain\)
Intermediate TetR dimer that has one specific site cleaved	\(cTetR:TetR\)
Intermediate TetR dimer that has one TetR cleaved binded to Cruzipain	\(cTetR:TetR:Cruzipain\)
Intermediate TetR dimer that has both specific sites cleaved	\(cTetR_{2}\)
Intermediate TetR that has its specific site cleaved	\(cTetR\)
TetR dimer binded to Cruzipain	\(TetR_2:Cruzipain\)
TetR binded to Cruzipain	\(TetR:Cruzipain\)

The following reactions were modelled:

$$ 1. DNA + RNAP \leftrightharpoons DNA:RNAP $$ $$ 2. DNA:RNAP \to DNA:RNAPa $$ $$ 3. DNA:RNAPa \to DNA + RNAP + mRNA $$ $$ 4. mRNA \to 0 $$ $$ 5. mRNA + Rib \to mRNA + Rib + Hirudin$$ $$ 6. Rib \to 0$$ $$ 7. DNA + TetR_2 \leftrightharpoons DNA:TetR_2 $$ $$ 8. DNA:TetR_2 \to DNA:TetRa_2 $$ $$ 9. TetR_2 \leftrightharpoons TetR + TetR$$ $$ 10. DNA:TetRa_2 + Cruzipain \leftrightharpoons DNA:TetRa_2:Cruzipain $$ $$ 11. DNA:TetRa_2:Cruzipain \to DNA + cTetR:TetR + Cruzipain $$ $$ 12. cTetR:TetR + Cruzipain \leftrightharpoons cTetR:TetR:Cruzipain $$ $$ 13. cTetR:TetR:Cruzipain \to cTetR_2 + Cruzipain $$ $$ 14. TetR_2 + Cruzipain \leftrightharpoons TetR_2:Cruzipain $$ $$ 15. TetR_2:Cruzipain \to cTetR:TetR + Cruzipain $$ $$ 16.TetR+Cruzipain \leftrightharpoons TetR:Cruzipain $$ $$ 17.TetR:Cruzipain \to cTetR + Cruzipain $$ $$ 18.cTetR:TetR \leftrightharpoons cTetR + TetR $$ $$19. cTetR_2 \leftrightharpoons cTetR + cTetR $$

These reactions were modelled in ODEs and simulated using the ode15s function at an absolute tolerance of 10-30 and a relative tolerance of 10-7.

1. Production of hirudin

For this simulation, we wanted to know how fast a cell-free system can produce hirudin. This is very important for us to determine if our system will be able to prevent blood coagulation. Blood clots in 5-10 minutes when taken out of the body, so our system will need to produce enough hirudin in less than that time.

To determine the production rate, we measured the time it takes to produce a threshold level of hirudin. We decided to use 1.3 μM of hirudin as the amount needed to prevent blood coagulation. (Markwardt, 1992)

By running our model with the following initial conditions:

• 50 nM of repressed DNA
• 50 nM of excess TetR dimer
• 360 pM of Cruzipain (Positive Test)
• 30 nM of RNAP
• 30 nM of Ribosomes

We only used repressed DNA in this simulation, because that will give us the rate of production of hirudin without considering the leakiness of the promoter. In other words, this should give us the slowest production rate in the scenario.

Results:

[Fig 1]

From the graph above, we see that it takes around 31 minutes to produce the threshold amount of hirudin. This delay is too long for the hirudin to effectively inhibit blood coagulation. In fact, running the blood coagulation model using this output revealed the presence of a thrombin spike at close to 8 minutes, indicating the initiation of blood coagulation. See the blood coagulation model for more details.

[Fig 2]

To see if we could increase the hirudin production rate, we tried changing the concentration of DNA.

[Fig 3]

We ran the simulations from 10 nM to 100 nM of DNA at 10 nM intervals. We see that increasing the amount of DNA can reduce the time taken to threshold, however, there is saturation at high concentrations of DNA. Further increasing the DNA concentration in our simulations revealed that it is unable to reduce the threshold time to lower than 20 minutes (at 3 μM of DNA), which is still too slow to inhibit blood coagulation.

Another initial condition that could be changed is the amount of excess TetR that is not bound to a promoter. To reduce chances of leakage, we assumed that the amount of TetR added to the kit will be higher than the amount of DNA. However, these excess TetR will competitively inhibit cruzipain cleavage of the TetR that is bound to the promoters. This can be seen in the graph below

[Fig 4]

[Fig 5]

We ran simulations at different concentrations of excess TetR from 0 nM to 50 nM, but found very little difference in the hirudin threshold time.

Finally, we ran a sensitivity analysis on the model parameters to determine the most suitable parameters to change to have the best impact on the hirudin production rate.

[Fig 6]

The vertical axis represents the different model parameters and the horizontal axis represents the relevant outputs of the model. The first five columns are the end-point (t = 3600s) concentrations of the respective species. The last three columns are metrics representing the dynamics of hirudin production. Initial hirudin Rate is taken to be the average rate over the first 10 seconds, which is useful in knowing the response delay. hirudin Rise Time is the time it takes to reach 90% of steady state (in this case, the end-point), which is a common way of measuring a system’s dynamic. The last is the just the maximum derivative of hirudin concentration.

The sensitivity of each parameter to each output is normalized to represent the ratio of the relative change in output to the relative change in the parameter. We see that changing the reaction kinetics of cruzipain will have a large effect on the final concentration of TetR, which is as expected. To find the best parameters to increase hirudin production rate, we see it is sensitive to the parameters k_TX, k_TL, ka_RNAP_DNA and K_Rib. Hence, changing the promoter strength or ribosome binding site (RBS) strength could increase the production rate.

However, a better method is to add amplification to our system.

2. Assumptions

Due to the slow rate of directly expressing hirudin. We decided to include an amplification step into our system. This is achieved by producing TEV instead, which will then cleave a sterically inactivated form of hirudin to release hirudin.

This is modelled by adding the following reactions and replacing the translation product with TEV.

[Fig 7]

The graph shows that using an amplification step will greatly increase the speed of the release of hirudin. The threshold time is reduced to 11 minutes from 31 minutes. Even though the production of TEV is slower as TEV has a longer sequence, it was still able to accelerate the release of hirudin.

As 11 minutes is very close to the common blood coagulation initiation time of 5-10 minutes, it is hard to tell if the blood will clot or not. Hence, we ran the blood coagulation model with this model’s output.

[Fig 8]

The lack of a distinct peak here indicates no initiation of blood coagulation, showing that an amplification step is effective.

Changing the Ribosome Biding Site Strength

To further increase the rate of production, we wanted to see what does changing the ribosome binding site will do to the production of Hirudin. We simulated the model with a Ribosome to RBS dissociation constant 2 orders of magnitude higher and lower than the current.

[Fig 9]

Lowering the Kd of the RBS increases its strength and increasing it decreases the strength of the RBS. However, we see that a weaker RBS has a larger effect on the production rate than a stronger RBS. Hence, our RBS seems to give very close to the maximum rate of hirudin production.

Stochastic Model – Extrinsic Noise

Most biochemical processes are stochastic in nature due to the uncertainty in the model’s parameters and the inherent unpredictability in the processes. Noise in biological systems can be classified into two groups – intrinsic noise and extrinsic noise. Intrinsic noise is the random fluctuations in reaction rates caused by the inherent stochasticity of biochemical reactions. It is more apparent in systems with a low concentration, where the low probability of collision creates large deviations in the rate of collisions. Due to our large reaction volumes of 30 μL, we have chosen to ignore the effects of intrinsic noise.

Extrinsic noise is different in that it is caused by the uncertainty or the variability in the model parameters. We modelled this by varying the parameters according to their measurement errors and running the deterministic model and plotting a histogram of the data.

References

Index	Reference
A	Dos Reis, F. C. G. et al. (2006) ‘The substrate specificity of cruzipain 2, a cysteine protease isoform from Trypanosoma cruzi’, FEMS Microbiology Letters, 259(2), pp. 215–220. doi: 10.1111/j.1574-6968.2006.00267.x.
B	Hillen, W. et al. (1983) ‘Control of expression of the Tn10-encoded tetracycline resistance genes. Equilibrium and kinetic investigation of the regulatory reactions’, Journal of Molecular Biology, 169(3), pp. 707–721. doi: 10.1016/S0022-2836(83)80166-1.
C	Karzbrun, E. et al. (2011) ‘Coarse-grained dynamics of protein synthesis in a cell-free system’, Physical Review Letters, 106(4), pp. 1–4. doi: 10.1103/PhysRevLett.106.048104.
D	Kleinschmidt, C. et al. (1988) ‘Dynamics of Repressor—Operator Recognition: The Tn10-encoded Tetracycline Resistance Control’, Biochemistry, 27(4), pp. 1094–1104. doi: 10.1021/bi00404a003.
E	Murray, D. N. A. M. G. and Biology, P. (2008) ‘Nucleic Acids Research’, Nucleic Acids Research, 36(19), pp. ii–ii. doi: 10.1093/nar/gkn907.
F	Stogbauer, T. R., 2012. Experiment and quantitative modeling of cell free gene expression dynamics, Munich: Ludwig–Maximilians–University.