Overview
Modelling is a really useful tool for understanding and predicting the result of the reactions between molecules. On this page is our perspective on modeling, using it in our project. It consists of the main aspects of design, mathematical and stochastic models of kinetic rates, but also includes the model of the binding affinity on our hardware. A more in depth-explanation of these models can be seen below!
Introduction
Gibbs free energy and stability
Figure 1. Gibbs free energy of a system where the reaction is spontaneous.
The change of DG is the value of energy that releases to the environment.
Gibbs free energy (G) of a system is the usable energy the system has. So DG is the change of Gibbs free energy and is the Gibbs free energy at a final state minus the Gibbs free energy at the initial state:
During a reaction, if DG of the system is negative, the reaction is spontaneous and will end in a product that is more stable than the reactants. The more negative DG is, the more stable the molecules are. In our case, the upper limit of free energy Gibbs for the molecules is -16 kcal/mol, in order to resist leakage in the absence of the initiator (Li et al. 2020).
Toehold mediated strand displacement
As a paying Slides customer you have access to the which allows you to add custom styles to your deck. By turning on the developer mode a new "class name" field will appear for any block that you focus.
Figure 2. Toehold mediated strand displacement reaction.
Taken from https://academic.oup.com/nar/article/41/22/10641/2437516#supplementary-data
The reactants of toehold mediated strand displacement are a double strand complex and a single strand, which is called invader(X). The complex consists of a strand called incumbent(Y) and the other one, the substrate(S), which includes a single-stranded overhang, the toehold. The invader strand is fully complementary, according to Watson-Crick rules, to the substrate, and these two strands can bind to each other using the toehold domain. This reaction is reversible as the toehold may dissociate. As the reaction proceeds, at the state of branch migration, the three-stranded complex of invader, incumbent and substrate, is moving back and forth. The reaction is completed when the incumbent is released. All in all, the displacement is thermodynamically driven forward by the net gain in base pairs due to the toehold. (Srinivas et al. 2013).
Molecules Design
In our project, we designed 2 different types of DNA hairpins. The tertiary structure of the first type of hairpin assimilates a Y structure molecule and it consists of a sequence complementary to miR-21, a sequence complementary to miR-10b and an initiator sequence. The release of the initiator sequence activates the second type of hairpins, the HCR hairpins. The HCR hairpins are DNA/RNA molecules. Our team designed 4 different HCR hairpins that, in the presence of the initiator and when thrown together, change secondary structures due to their complementarity.
Our goal was to design stable molecules that do not unfold in the absence of the initiator. The stability of the molecules is determined through the Gibbs free energy that characterizes them. If and only if the two microRNAs are present, the Y molecule destabilizes and the initiator is released. The initiator is a single strand of DNA, which commences the process of Hybridization Chain Reaction amplification (HCR).
Structure
The design of our molecules is based on complementary sequences of the inputs in every step. In theory, any guanine can bind to any cytosine, and any adenine can bind to any thymine, but in fact the reverse reaction is so fast, that the expected lifetime of an isolated one-base-pair-binding is very short. A stretch of several consecutive nucleotides must be complementary for the bound state to exist as an intermediate for further reaction at longer time scales.
Y-structure Hairpins
In the Y structure, the Y1 and Y3 sequences are complementary to mir 21 and mir 10b, respectively. In strand Y1 we added some unpaired bases, which is called bulge (Ma et al. 2019).
Figure 3. Structure of Y-hairpin
According to Leontis et al. (1991), the hairpin’s stability increases by adding unpaired bases, the bulge, on three-way DNA junctions, as in presence of metal ions, the structure of the three-way junctions bend into a well-defined structure.
Based on the sequences of the hairpins, the folded structures appear in two stereochemically different isomers. Most of the sequences that were tested adopt the isomer I conformation. The bulged DNA molecules are ion- and sequence-dependent. (Welch et al. 1995)
Figure 4. Isomers of Y-hairpins. Taken by Welch et al. 1995
HCR hairpins
Figure 5. HCR hairpins(* domains are complementary to domains with the same letter, e.g a and a*)
As seen above, HCR hairpins can be separated in three main parts: loop, stem and toehold. In our project, the HCR hairpins bind to each other because of their complementarity in parts.
- The initiator, which, as we explained earlier, begins the Hybridization Chain Reaction is binding to H1 hairpin, as it is complementary to its toehold and its stem(a* and b*).
- H1 hairpin’s loop and stem part(c and b), which are as a straight strand of DNA, are binding to H2 hairpin’s toehold and stem(c* and b*), respectively.
- H3 and H4 hairpins bind to H1-H2-initiator structure, in the same way as H1 and H2 were bound.
- As the HCR hairpins bind to each other, therapeutic molecules siRNAs are released, in order to form siRNA duplexes.
Design Parameters:
- Toehold length: The main type of reaction between molecules in our project is toehold mediated strand displacement. In more detail, this reaction is initiated by a short single-stranded overhang region, known as a toehold, which can greatly increase the rate of strand exchange. In toehold-mediated strand displacement reactions, there is a strong coupling between the kinetics and the thermodynamics. Thus, to accelerate the strand displacement reaction, the invading toehold must be made stronger thermodynamically by increasing the toehold length.
- Loop length: It is known that hairpins with a loop with bigger length and little stem are more stable than little loop and bigger stem.
- Illegal sequences of iGEM: We have tested through our software on design algorithms some sequences that must not be included on our haiprins, because they act as recognition sites on some enzymes.
HCR step Stability
HCR hairpins are hybrid molecules, as they consist of DNA and RNA. DNA/RNA interactions,as they observed between the HCR hairpin and siRNA or antisense sequence, are known to be less stable than DNA/DNA and RNA/RNA, as well.
DG for hybrid molecules, like HCR hairpins, is calculating by the equation:
So, DG initial on HCR step is:
As for the final state at HCR, Gibbs free energy is:
where: 
The difference at DG from the final to initial state is:
DDG is negative, as 
and the product from HCR reaction is more stable than the reactants. All of these equations and Gibbs free energy are calculated through design software on github, which consists of nupack package (Gong et al. 2021).
Calculation of Gibbs free energy of our parts
We have used the code that is available on github, where we have incorporated the last updated nupack package, to determine the Gibbs free energy of the reactants and the product of every reaction. We have set the temperature at 37 o Celsius.
Table 1. Free energy Gibbs of all the reactants molecules
All the results can be found in the supplementary pdf on the end of this page.
Intuitive Energy Landscape Model
Modeling main aspects of the project will be helpful, not only to understand the parameters of the molecular reactions, but also to predict whether these reactions will take place. Modeling the toehold mediated strand displacement reaction bases on the intuitive energy landscape (IEL) model. It is an analysis on the biophysical basis of the molecular mechanism of toehold-mediated strand displacement, which regulates the control of kinetics of molecular rearrangements.
Models
Gibbs free energy model
IEL models the free energy Gibbs of a virtual box of volume V, which consists of the two invalid molecules, the reactants, in a concentration u per volume V. The two parameters that characterize the IEL model are DGs ,“the effective sawtooth amplitude”and DGp ,the plateau height.
- The first parameter, DGs, is referring to the state where toehold is fully bound and branch migration transpires. The adjustability of its unit can control if branch migration and toehold melting rates act independently.
- The plateau height, DGp, the second parameter, captures how the free energy of branch migration intermediates could vary with the structure of the branch migration junction and depends on entropic or electrostatic effects. In the IEL model, a restriction for the plateau height is that it can’t be equal to zero, in order to measure hybridization, fraying and branch migration rates, and hence strand displacement rates, to be matched.
Figure 6.Gibbs free energy box of the reaction of toehold mediated strand displacement according to the IEL model.
This figure refers to the reaction of Y-structure with the miRNA-21. The different colours refer to the different concentration the IEL model has tested.
Table 2. Description of the different states of Gibbs free energy box figure.
The free energy for the first step (state A to B), the initial binding, consists of the energy for the binding at a standard concentration, u0=1M , DGassoc, and the energy for the correction to the concentration, DGvolume.
So, the equation for the initial binding is:
DGinitial =DGassoc+ DGvolume,
where:
DGvolume=R*T*ln(u0/u)
Assumptions:
- It is a sequence independent model, but length dependent
- The reactants are referring as unfolded structures
- The units for the parameters are using by the results according to the paper
Kinetic rates model
We have used the figures of free energy to understand the steps on the toehold mediated strand displacement reaction. The transition rates, as they satisfy the thermodynamic Boltzmann equilibrium, are calculating by the equation:
,where kij is the transition rate from state i to j
The kinetic rates, which are independent from the free energy of the system, are the unimolecular kinetic rate and the bimolecular rate, for unimolecular and bimolecular transitions, respectively. Unimolecular transitions are referring to all states, except A to B and D to E, which can be characterized as bimolecular transitions.
For unimolecular transitions: kij = kuni
and thus,
On the other side, for bimolecular transitions, the join step of the compounds, i to j state, with the kinetic rate:
and thus,
Based on the hybridization rate constant fitted by Zhang and Winfree, which is referring to strands with long toehold, as there are in the project, kbi=3*106/M/s .
The unimolecular reaction rate is calculating by the equation, using Metropolis Method:
Τhe overall effective rate constant, keff, is given by the equation, where h is the toehold length and bis the branch migration length, is:
where n, is given by the equation:
The average amount of time needed for branch migration step in either direction is called the average branch migration step time (tbm) and is calculated by the equation:
where kbm, the average rate constant of the branch migration step is:
Probabilities between the transitions
Really important to understand the different states of the toehold mediated strand displacement reaction is the estimation of the probabilities between these transitions.
So, the probability from state i to i+1 is:
and the probability from state i to i-1:
The probability of fully bound toehold is calculating by the equation:
where,
The probability of the displacement of the incumbent before the dissociation of the toehold, which is fully zipped up, is:
(Srinivas et al. 2013).
Results from using IEL on our project
We tested the IEL model on our parts, using the toehold and branch migration length of their sequences. The figures about the Gibbs free energy and the kinetic rates and the table of the probabilities are below.
For example, for the first step, the Y-structure hairpin reacts with the miRNA-21 and the miRNA-10b and produces the sequence of the Initiator. In the first reaction with miRNA-21, the toehold length is 8 and in the second reaction with miRNA-10b is 6. The branch migration length is 14 and 17 in the first and second reactions, respectively. Running the code of IEL for these reactions, we got the results below.
Running every step for different concentrations, to determine the one that leads to better results. So, for the reaction of the Y-structure hairpin with miRNA-21, for concentration 0.1-4.6 M with a step of 0.5 M we got the figures, for the DG box(figure 6), the kforward, the kback, the table of pforward and pback
Figure 7. The rate constants of the forward and the backward transitions of every step
Table 2. The probabilities of the forward and the backward initial step according to different concentrations.
As we can see from the results above the rate constant for the forward reactions for high concentration is bigger than the rate constant for the backward reaction, although it is increasing by the number of steps on toehold reaction. As for the probabilities values are really high when the concentration is low and as it increased the probability on the forward initial step is much higher than for the probability to the backward reaction.
All the results can be found in the supplementary pdf on the end of this page.
Uncertainty Analysis
The uncertainty analysis aims to describe the uncertainty of variables that are quantified by a model. The propagation of uncertainty is referring to a system of variable’s uncertainty, which investigates statistical calculation of the uncertainty of setted variables to aim for an accurate value of uncertainty. There are multiple methods to carry out this analysis (Saltelli et al. 2008).
Figure 8. Uncertainty analysis methodology, taken by Coleman et al. 2009
We chose Monte carlo method simulation, which basic methodology, as it is shown in the figure above, is:
i. Input the value, the elements of error source, the probability distribution function and the range of uncertainty related to each variable. According to the central limit theorem, even if the original random variables themselves are not normally distributed, when the independent random variables are added together, their correctly normalized total tends toward a normal distribution
ii. Define the number of trials, as they are the number of iterations on monte carlo simulation
iii. For every iteration, random value of each ε is selected and Xi, simulated measured value
iv. Calculate the standard deviation of Xi
v. Check for convergence (plot sMCM)
vi. Calculate the expanded uncertainty (Coleman et al. 2009).
As we have determined the effective rate, keff, of each transition on the Y-structure hairpin and the HCR hairpins, we run the code of IEL model, we define for every one the value of uncertainty, setting the uncertainty of the independent variables as 5%, (at a 95% level of confidence) . The figures below refer to the first transition of Y hairpin with the miRNA 21.
Figure 9. Results of the uncertainty analysis on the initial transition step of Y-structure with the miRNA-21.
In the figures below, it is shown the distribution of the result keff at a probabilistically symmetric coverage interval, the best value of keff considering the uncertainty of all the independent variables and the standard deviation smcm.
All the results can be found on the supplementary on the end of this page.
Gillespie Stochastic Simulation
As we have calculated with the uncertainty analysis the value of uncertainty of kinetic rates, we use a stochastic model in order to simulate the reaction and how the concentration of the molecules fluctuates. Using stochastic simulation models, we consider that collisions in a system of reactions of molecules occur in an random manner. The reaction for Gillespie Algorithm is shown below.
The Gillespie Stochastic algorithm, which is used to simulate chemical or biochemical systems of reactions, creates a statistically correct trajectory of a stochastic equation system for which the reaction rates are known. It is based on a random-walk process, with a parameter of randomness between 0 and 1. (Nezar et al. 2022).
On Gillespie stochastic algorithm, according to the reactions above, we firstly state the initial variables of the concentration of molecules, the rate constants and fill the matrix, according to stoichiometries of the reaction. Also, we input the simulation time range. So using the Gillespie model, we estimate the concentration of the molecules of our project as they are shown on the figures below.
Figure 10. Results of Y-structure step reaction
Figure 11. Results of HCR reaction molecules
So, as we saw on figure the concentration of Y2-Y3 and the initiator as they produced, they act simultaneously as reactants of the next transition. The molecules of the initiator have been used in the step of HCR. In the figure, according to the rate constants of the transitions of HCR step, the molecules of H1-Initiator as it produces it reacts concurrently in the next reaction. Also, the H3-H2-H1 initiator, because of the value of rate constants, has a high concentration, like the HCR product, which is the final product.
Hardware
Modeling hardware parts is something that we came up with as we were looking for its design. It helped us understand more its characteristics and guide us on the densities that we used in the experiments.
Our hardware is an electrochemical DNA sensor using gold microelectrodes, on which are immobilized specific sequences of DNA. These sequences are complementary to the target strands that we want to detect. Thus, it includes two main states, one with free and one with target-bound receptors. The relative fraction of target-bound receptors can be described by the Langmuir adsorption isotherm:
where θ is referring to the target-bound fraction, T is the concentration of the target, the probe that is immobilized on the electrochemical sensor, and KD the dissociation constant of the receptor.
The assumptions that have been used in this model are (Pellitero et al. 2019):
- only one binding site is available at each receptor
- all probes on the surface have the same binding affinity for the target
- the DNA probes do not interact with each other
- the extent of target binding is limited to a single monolayer
In the development of DNA sensors, a major restriction is the ability of immobilization of DNA probes to gold surface, optimizing its density.The signal of the sensor is getting higher as the value of target DNA, which is hybridized with probes DNA, is increasing. According to literature research, the limit of density of probes where the signal is maximum is 4 × 1012 probes/cm2. At higher densities, the value of hybridized targets decreases. Hybridization kinetics are referred to as fast when the amount of probes is 3 × 1012 probes/cm2 or below. Also, the density of DNA on the surface is known to be sequence and length dependent.
Another parameter that affects the signal of the DNA sensor is the geometric characteristics of the surface of the microelectrodes. This connection is described by the value of roughness factor, which can be equal from 1 to 3, according to different types of gold microelectrodes and the use of different protocols (Keighley et al. 2008).
Thermistor Resistance
One of our hardware elements is the temperature control system, a temperature sensor called a thermistor. A thermistor is a semiconductor, in which when the temperature rises, the resistance of the thermistor falls and the resistance rises as temperature falls (NTC thermistor model. 2022 Oct. 12).
The resistance of the thermistor can be calculated by the equation:
R = Ro exp( Beta/T - Beta/To)
where:
- R: the resistance of the thermistor at T(K)
- Ro: the normal resistance of the thermistor at To(K)
- Beta: the thermistor material constant
- To: the temperature where T has been measured
By the model of the thermistor we understand its operation; thus we incorporate these variables on our Arduino code.
Other MDS
Collaboration with iGEM Mohali
As our model is sequence independent, in the context of our collaboration with iGEM Mohali, they helped as they have run on their molecular dynamics simulator our sequences of the Y-structure hairpin and they send us the results, which are shown in the figure below.
Figure 12. Stability of the complex of Y hairpin with the microRNAs according temperature
Their Molecular Simulation of RBD, which is a sequence dependent model, simulates the binding to a functional aptamer. From this figure, we understand that this complex is stable at a temperature of 300 K.
References
Coleman H.W. and Steele W.G. (2009) “Experimentation, Validation, and Uncertainty Analysis for Engineers, Third Edition”, John Wiley & Sons, Inc. ISBN: 978-0-470-16888- 2
Pellitero, M. A., Shaver, A., & Arroyo-Currás, N. (2019). Critical review—approaches for the electrochemical interrogation of DNA-based sensors: A critical review. Journal of The Electrochemical Society, 167(3), 037529. https://doi.org/10.1149/2.0292003jes
Gong, X., Wang, H., Li, R., Tan, K., Wei, J., Wang, J., Hong, C., Shang, J., Liu, X., Liu, J., & Wang, F. (2021). A smart multiantenna gene theranostic system based on the programmed assembly of hypoxia-related siRNAs. Nature communications, 12(1), 3953. https://doi.org/10.1038/s41467-021-24191-9
Keighley, S. D., Li, P., Estrela, P., & Migliorato, P. (2008). Optimization of DNA immobilization on gold electrodes for label-free detection by electrochemical impedance spectroscopy. Biosensors & bioelectronics, 23(8), 1291–1297. https://doi.org/10.1016/j.bios.2007.11.012
Li, S., Li, P., Ge, M., Wang, H., Cheng, Y., Li, G., Huang, Q., He, H., Cao, C., Lin, D., & Yang, L. (2020). Elucidation of leak-resistance DNA hybridization chain reaction with universality and extensibility. Nucleic acids research, 48(5), 2220–2231. https://doi.org/10.1093/nar/gkaa016
Leontis, N. B., Kwok, W., & Newman, J. S. (1991). Stability and structure of three-way DNA junctions containing unpaired nucleotides. Nucleic acids research, 19(4), 759–766. https://doi.org/10.1093/nar/19.4.759
Ma, X., Chen, X., Tang, Y., Yan, R., & Miao, P. (2019). Triple-Input Molecular AND Logic Gates for Sensitive Detection of Multiple miRNAs. ACS applied materials & interfaces, 11(44), 41157–41164. https://doi.org/10.1021/acsami.9b16812
Nezar (2022). Gillespie Stochastic Simulation Algorithm (https://github.com/nvictus/Gillespie), GitHub. Retrieved October 11, 2022.
Saltelli A., Ratto M., Andres T., Campolongo F., Cariboni J., Gatelli D., Saisana M. and Tarantola S. (2008) “Global Sensitivity Analysis. The Primer”, John Wiley & Sons, Ltd. ISBN: 978-0-470-05997-5
Srinivas, N., Ouldridge, T. E., Sulc, P., Schaeffer, J. M., Yurke, B., Louis, A. A., Doye, J. P., & Winfree, E. (2013). On the biophysics and kinetics of toehold-mediated DNA strand displacement. Nucleic acids research, 41(22), 10641–10658. https://doi.org/10.1093/nar/gkt801
Welch, J. B., Walter, F., & Lilley, D. M. (1995). Two inequivalent folding isomers of the three-way DNA junction with unpaired bases: sequence-dependence of the folded conformation. Journal of molecular biology, 251(4), 507–519. https://doi.org/10.1006/jmbi.1995.0452
