ResiSense
Finding resistant pathogens one species at a time
🔍
Use ↑, ↓, and Enter to navigate the results.
Querying took 0ms. Redrawing took 0ms.
Modelling
Modelling
Modelling

This page details our dry lab work on the G-quad-haemin complex. To download our various resources, see the resources page.

Introduction

In our project, G-quadruplex DNAzyme (G-quad) is the last step in our detection mechanism, which produces a visible colour change when the target gene is detected. The colour change is due to the peroxidase-mimicking activity when the G-quadruplex DNAzyme is bound with haemin (G-quad-haemin), which catalyses the oxidation of TMB (colourless) to TMB⁺ (blue) in the presence of H₂O₂, as shown in Equation 1 [1]. The literature has shown that the binding of haemin to G-quad is crucial to the catalytic activity of G-quad [2], and G-quadruplex is more stable in a mixture of water and dimethyl sulphoxide (DMSO) [3].

Since two of the reagents in the reaction, haemin and TMB, also have better solubility in DMSO than in water [4], [5], it is postulated that adding more DMSO to the reaction mixture may help G-quad work faster. To investigate this, a molecular dynamics (MD) approach has been developed to calculate the binding energy of haemin to G-quad at different concentrations of DMSO in water.

Equation 1: Oxidation of TMB into TMB⁺ by G-quad-haemin complex

Equation 1: Oxidation of TMB into TMB⁺ by G-quad-haemin complex

Methodology

To obtain the free energy of binding of haemin to G-quad, alchemical transformation of the haemin is performed, and the free energy of binding is calculated according to the thermodynamic cycle, as shown in Figure 1. The software being used is GROMACS [6].

Figure 1: Thermodynamic cycle for calculating the binding energy of haemin to G-quad

Figure 1: Thermodynamic cycle for calculating the binding energy of haemin to G-quad

Figure 1: Thermodynamic cycle for calculating the binding energy of haemin to G-quad

Topology Preparation

The topology is prepared using the AlphaFold 3 model [7], with one DNA input (5'-GGGTGGGTGGGTGGGTC-3') and one ligand input (HEM - Haem). It is noted that haemin and haem are different, with haemin having an extra chloride ion coordinated with the iron ion (Figure 2). Since haemin and haem have very similar structures, we believe that haem is a sufficient approximation to haemin. The model generated by AlphaFold 3 with the highest confidence is taken as the topology for simulation.

Figure 2: Structure of haemin (left) and structure of haem (right) Figure 2: Structure of haemin (left) and structure of haem (right)

Figure 2: Structure of haemin (left) and structure of haem (right)

Topology Verification

We have visually inspected the model generated to ensure that AlphaFold 3 provided a reasonable starting topology for simulation. It was observed that the generated model follows the end-stacking binding mode reported in the literature [8], as shown in Figure 3 and Figure 4.

Figure 3: Folded G-quadruplex DNAzyme generated by AlphaFold 3

Figure 3: Folded G-quadruplex DNAzyme generated by AlphaFold 3
Note that the G-quadruplex structure in the middle of the DNAzyme could be clearly seen.

Figure 4: Folded G-quadruplex DNAzyme bound with haem

Figure 4: Folded G-quadruplex DNAzyme bound with haem
Note that the haem stacks on top of the G-quadruplex, which is consistent with the end-stacking binding mode.

Choice of Force Field

The amber14sb_OL15 force field [9], [10] is chosen as the force field to be used in the simulation, as this force field has been found to better resemble experimental measurements than other force fields, which would provide a more realistic simulation [11].

Since DMSO and haem are not "native" molecular species to the force field, they have to be parametrised before the simulation could proceed.

DMSO Parametrisation

The DMSO molecule was parametrised by the ACPYPE server [12], [13]. The topology of DMSO was submitted to the server as a .pdb file, and the server produces the .itp and .gro file required for building the topology.

Haem Parametrisation

Since the center of haem contains an Fe³⁺ ion, haem cannot be parametrised by the ACPYPE server as the Generalised Amber Force Field (GAFF) does not contain the Fe³⁺ ion [14]. To parametrise haem in a way that is compatible with amber14sb_OL15 force field, the Metal Center Parameter Builder (MCPB.py) is used for preparing the topology [15].

First, the haem molecule is passed to MCPB.py to obtain input files for GAMESS-US [16]. Then, quantum calculations were done using GAMESS-US in order to optimise the geometry and compute the charge density. The wave function is computed using B3LYP/6-31G level of theory*, which is sufficient for haem [17]. Restrained Electrostatic Potential (RESP) [18] was applied to compute the partial charges on different atoms, and the Seminario method [19] was used to compute the bonded parameters. Lastly, the .frcmod file were converted to .itp and .gro using the LEAP program and ACPYPE.

Topology Assembly

While DMSO and haem were parametrised in previous steps, they are not a part of the residue database. Therefore, they have to be handled separately from the G-quad.

First, the .pdb file generated by AlphaFold 3 was separated into two .pdb files, one for G-quad and one for haem. The G-quad is converted into GROMACS topology using the pdb2gmx command. For the haem, the editconf command is used to convert .pdb to .gro. The .gro coordinates were manually copied to the .gro file of the G-quad .gro file. The .top file was accordingly modified to include the topology for haem and the restraints for it.

After that, if DMSO is used as a co-solvent, the insert-molecules command is invoked to add DMSO molecules into the topology. solvate is used to add water to the system. The SPC/E water model [20] is used. To neutralise charges, genion is used to add ions to the topology. K⁺ ions are added to balance negative charges, and Cl⁻ ions are added to balance positive charges, as K⁺ ions are present in the experimental setup and were also shown to have the best effect on stabilising G-quad structures [21].

System Equilibration

To arrive at a reasonable starting structure and to have correct sampling during the production MD, the system must be well-equilibrated before the production MD runs. Energy minimisation was performed on the solvated structure using the steepest gradient descent method [22]. NVT equilibration was performed on the energy-minimised structure. No additional thermostats were necessary since the Langevin dynamics integrator was used for the simulation, which could also regulate the system temperature. The system temperature was set at 298K. NPT equilibration was performed on the NVT ensemble, using the Parrinello-Rahman barostat [23] to regulate the system pressure. The pressure is set to be 1 bar. The NVT and NPT equilibration were run for 100ps. The potential, temperature, pressure, and density over time were plotted, and convergence was observed in the graphs.

Production MD

A total of 41 production MDs were run for each solvent configuration, in which the charge interaction of haem is gradually decoupled first, and then the Lennard-Jones potentials are also gradually decoupled. This allows us to obtain the ∂H/∂λ curve, which can then be used to obtain ΔG using the Bennett Acceptance Ratio (BAR) method [24].

The charge lambdas are first decoupled in steps of 0.05 per iteration, and after the charge interactions are decoupled, the Lennard-Jones potentials are decoupled in steps of 0.05 per iteration. To ensure sufficient ensemble sampling at the specified lambda values, 2.5ns of the simulation was run for each combination of lambda, which is around two times what is deemed adequate [25]. The free energy was obtained from ∂H/∂λ using bar.

Calculation of ΔG of G-quad-haem Bond

To calculate the ΔG of haem binding to G-quad, one may refer to the thermodynamic cycle outlined in Figure 1. Two separate simulations were set up, one with G-quad bounded with haem and another with haem only. Two systems were solvated in the same DMSO/water mixture concentration, and haem was alchemically transformed through the process outlined above.  

Results

The simulation results are summarised below. ΔG are given in kJ/mol and the +/- values represent the standard deviation of ΔG.

DMSO concentration (v/v) % ΔG in G-quad-haem ΔG in haem only ΔG binding
0 2598.74 +/- 15.45 2633.05 +/- 5.10 34.41 +/- 16.27
5 2577.23 +/- 3.00 2592.05 +/- 7.48 14.82 +/- 8.06
10 2577.19 +/- 11.07 2615.68 +/- 2.05 38.49 +/- 11.26
15 2562.85 +/- 4.34 2604.41 +/- 2.51 41.56 +/- 5.01
20 2510.37 +/- 17.08 2626.85 +/- 1.66 116.48 +/- 17.16
25 2543.50 +/- 11.67 2568.09 +/- 2.83 24.59 +/- 12.01
30 2482.45 +/- 9.99 2600.48 +/- 3.56 118.03 +/- 10.61
35 2492.22 +/- 13.88 2594.15 +/- 5.82 101.93 +/- 15.05
40 2564.71 +/- 22.22 2586.21 +/- 3.90 21.50 +/- 22.56
45 2479.50 +/- 18.14 2558.66 +/- 3.21 79.16 +/- 18.42
50 2509.47 +/- 15.19 2524.21 +/- 9.77 14.74 +/- 18.06

To visualise the results, the free energies are plotted against DMSO concentration (v/v) %, as shown in Figure 5 and Figure 6.

Figure 5: Change in free energy (ΔG) of decoupling of haem in G-quad-haem and haem only

Figure 5: Change in free energy (ΔG) of decoupling of haem in G-quad-haem and haem only

Figure 6: Change in free energy (ΔG) of binding of haem in G-quad-haem

Figure 6: Change in free energy (ΔG) of binding of haem in G-quad-haem

Discussions

From Figure 5, we could observe that the free energy change of decoupling of haem only in solvent decreases as DMSO concentration increases. This is in agreement with the literature, where the solubility of haemin increases as DMSO concentration increases. In Figure 6, we observed a general trend of increasing binding energy of haem to G-quad, indicating a prediction of decreasing binding affinity to G-quad as the DMSO concentration increases.

Some data points do not follow the trend: the binding energies at 25% and 40% DMSO concentration (v/v). We believe this is due to the change in G-quad conformation during the simulation, which skewed the free energy values. Through the inspection of RMSD value over time of simulations (see Supporting Information) and comparing with the non-perturbed simulation (λ = 0), we conclude that the G-quadruplex DNAzyme has changed conformation in some of the simulations, leading to inaccurate free energy values. Potential solutions are discussed below.

We have also noticed that the standard deviations of the change in free energy for G-quad-haem decoupling are large. Upon inspection of the gmx bar output, we noticed that the relative entropies between some of the lambdas were significant, indicating inadequate phase space sampling. This suggests the need for more sampling in those intervals.

Limitations and Uncertainties

While the simulations have provided valuable insights in terms of the affinity of haem to G-quad, we are also aware of the following limitations and uncertainties.

Inability to Simulate Peroxidase Activity of G-quad-haem in Various Solvent Mixtures

While molecular dynamics simulations may provide information about the binding affinity between G-quad and haem, they cannot provide information on the peroxidase activity of G-quad-haemin due to their inability to simulate bond breaking and forming, which are the basis of chemical reactions. Such processes would lead to topology modifications, which is impossible in pure molecular mechanics simulations. Also, in molecular mechanics simulations, quantum effects on the electrons are neglected, and molecules are assumed to be at the ground state. In contrast, most reactions require electrons to be excited or exchange electrons between different species. Resolving such issues would demand the model to take quantum effects into account. We suggest that quantum mechanics/ molecular mechanics (QM/MM) hybrid simulations may be able to fill this gap.  

Unrestrained G-quad Skews Free Energy Values

Upon inspection of the G-quad-haem structure generated by the AlphaFold 3 model, we noticed that there are hand-like structures formed by the last two nucleotides of the G-quad sequence (TC-3') that seem to hold onto the haem ligand and keep it in place. However, as the interaction between haem and the surroundings was decoupled, the interactions between haem and the hand region weakened, which may cause the hand region to wiggle freely. This increases the number of possible microstates in the system, increasing the entropy and leading to skewed free energy values. This may be solved by "confine and release" free energy simulations, which have been proposed in literature [26].

Binding Affinity may not be a Universal Indicator to G-quad-haemin Peroxidase Activity

While there was literature suggesting that the binding affinity of haemin to G-quad may be a predictor for the peroxidase activity of the G-quad-haemin DNAzyme [8], we have also found literature supporting the opposing view [27]. Hence, the hypothesis that the free energy of binding between G-quad and haemin is a predictor for peroxidase activity may not be universally true.

G-quad-haem Activity may not be the Most Crucial Factor Affecting Reaction Rate

Although we focused on evaluating the stability of the G-quad-haem complex in different solvent mixtures, one must remember that the ultimate goal for optimising solvent mixture is to maximise the reaction rates, such that the test kit develops a more evident colour change more rapidly. The stability of the G-quad-haem complex only evaluated the catalyst in the reaction, and the effects of solvents on the reactants' stability may need to be accounted for in the simulations. The solvent may increase the reaction rate by affecting the stability of intermediate and final products and vice versa [28]. It is suggested that the whole reaction, including its intermediary products, may be qualitatively evaluated in terms of their stability before confirming whether the increased or decreased stability of G-quad-haem matters for the reaction.

Conclusion

Despite the presence of some inaccuracies due to the factors outlined above, our preliminary modelling results indicate a promising potential that the binding energy of haem to G-quad is the lowest in fully aqueous solutions. This exciting discovery could pave the way for a new method to predict G-quad-haem complex catalytic capability in different solutions, a possibility that can be further explored in the lab.

To help future iGEM teams in modelling molecular systems involving the haem residue, we have uploaded the parametrised haem topology for other teams to use.

References

  1. Z. Zhang, W. Zhao, C. Hu, Y. Cao, Y. Liu, and Q. Liu, "A Convenient and Label-Free Colourimetric Detection for L-Histidine Based on Inhibition of Oxidation of 3,3′,5,5′-Tetramethylbenzidine-H₂O₂ System Triggered by Copper Ions", vol. 9, p. 773519, Nov. 2021, doi: 10.3389/fchem.2021.773519.
  2. Y. Cao et al., "Investigation and improvement of catalytic activity of G-quadruplex/haemin DNAzymes using designed terminal G-tetrads with deoxyadenosine caps", vol. 11, no. 26, pp. 6896–696, Jul. 2020, doi: 10.1039/d0sc01905d.
  3. N. Tariq, T. Kume, L. Luo, Z. Cai, S. Dong, and R. B. Macgregor, "Dimethyl sulphoxide (DMSO) is a stabilising co-solvent for G-quadruplex DNA", vol. 282, p. 106741, Mar. 2022, doi: 10.1016/j.bpc.2021.106741.
  4. R. Stiebler, A. N. Hoang, T. J. Egan, D. W. Wright, and M. F. Oliveira, "Increase on the Initial Soluble Haem Levels in Acidic Conditions Is an Important Mechanism for Spontaneous Haem Crystallisation In Vitro", vol. 5, no. 9, p. e12694, Sep. 2010, doi: 10.1371/journal.pone.0012694.
  5. A. Frey, B. Meckelein, D. Externest, and M. A. Schmidt, "A stable and highly sensitive 3,3′,5,5′-tetramethylbenzidine-based substrate reagent for enzyme-linked immunosorbent assays", vol. 233, no. 1, pp. 47–56, Jan. 2000, doi: 10.1016/S0022-1759(99)00166-0.
  6. M. Abraham et al., "GROMACS 2024.3 Source code". Zenodo, 29 Aug, 2024. doi: 10.5281/zenodo.13456374.
  7. J. Abramson et al., "Accurate structure prediction of biomolecular interactions with AlphaFold 3", vol. 630, no. 8016, pp. 493–500, Jun. 2024, doi: 10.1038/s41586-024-07487-w.
  8. N. Alizadeh, A. Salimi, and R. Hallaj, Haemin/G-Quadruplex Horseradish Peroxidase-Mimicking DNAzyme: Principle and Biosensing Application, vol. 170. Switzerland: Springer International Publishing AG, 2017, pp. 85–106.
  9. M. Zgarbová, J. Šponer, M. Otyepka, T. E. Cheatham, R. Galindo-Murillo, and P. Jurečka, "Refinement of the Sugar–Phosphate Backbone Torsion Beta for AMBER Force Fields Improves the Description of Z- and B‑DNA", vol. 11, no. 12, pp. 5723–5736, Dec. 2015, doi: 10.1021/acs.jctc.5b00716.
  10. M. Zgarbová et al., "Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles", vol. 7, no. 9, pp. 2886–2902, Sep. 2011, doi: 10.1021/ct200162x.
  11. Li N., Gao Y., Qiu F., and Zhu T., "Benchmark Force Fields for the Molecular Dynamic Simulation of G-Quadruplexes", vol. 26, no. 17, pp. 5379-, Sep. 2021, doi: 10.3390/molecules26175379.
  12. Sousa da Silva A. W. and Vranken W. F., "ACPYPE - AnteChamber PYthon Parser interfacE", vol. 5, no. 1, pp. 367--367, Jul. 2012, doi: 10.1186/1756-0500-5-367.
  13. L. Kagami, A. Wilter, A. Diaz, and W. Vranken, "The ACPYPE web server for small-molecule MD topology generation", vol. 39, no. 6, Jun. 2023, doi: 10.1093/bioinformatics/btad350.
  14. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, "Development and testing of a general amber force field", vol. 25, no. 9, pp. 1157–1174, Jul. 2004, doi: 10.1002/jcc.20035.
  15. P. Li and K. M. Merz, "MCPB.py: A Python Based Metal Center Parameter Builder", vol. 56, no. 4, pp. 599–604, Apr. 2016, doi: 10.1021/acs.jcim.5b00674.
  16. Barca G. M. J. et al., "Recent developments in the general atomic and molecular electronic structure system", vol. 152, no. 15, p. 154102, Apr. 2020, doi: 10.1063/5.0005188.
  17. Autenrieth F., Tajkhorshid E., Baudry J., and Luthey‐Schulten Z., "Classical force field parameters for the haem prosthetic group of cytochrome c", vol. 25, no. 13, pp. 1613–1622, Oct. 2004, doi: 10.1002/jcc.20079.
  18. C. I. Bayly, P. Cieplak, W. Cornell, and P. A. Kollman, "A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model", vol. 97, no. 40, pp. 10269–10280, Oct. 1993, doi: 10.1021/j100142a004.
  19. J. M. Seminario, "Calculation of intramolecular force fields from second-derivative tensors", vol. 60, no. 7, pp. 1271–1277, 1996, doi: 10.1002/(SICI)1097-461X(1996)60:7<1271::AID-QUA8>3.0.CO;2-W.
  20. H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, "The missing term in effective pair potentials", vol. 91, no. 24, pp. 6269–6271, Nov. 1987, doi: 10.1021/j100308a038.
  21. S. Haider, "Computational Methods to Study G-Quadruplex–Ligand Complexes", vol. 98, no. 3, pp. 325–339, Sep. 2018, doi: 10.1007/s41745-018-0083-3.
  22. E. J. Haug, J. S. Arora, and K. Matsui, "A steepest-descent method for optimisation of mechanical systems", vol. 19, no. 3, pp. 401–424, Jul. 1976, doi: 10.1007/BF00941484.
  23. Parrinello M. and Rahman A., "Polymorphic transitions in single crystals: A new molecular dynamics method", vol. 52, no. 12, pp. 7182–7190, Dec. 1981, doi: 10.1063/1.328693.
  24. C. H. Bennett, "Efficient estimation of free energy differences from Monte Carlo data", vol. 22, no. 2, pp. 245–268, Oct. 1976, doi: 10.1016/0021-9991(76)90078-4.
  25. R. S. Rathore, P. Aparoy, P. Reddanna, A. K. Kondapi, and M. Rami Reddy, "Minimum MD simulation length required to achieve reliable results in free energy perturbation calculations: Case study of relative binding free energies of fructose-1,6-bisphosphatase inhibitors", vol. 32, no. 10, pp. 2097–2103, Jul. 2011, doi: 10.1002/jcc.21791.
  26. D. L. Mobley, J. D. Chodera, and K. A. Dill, "Confine-and-Release Method: Obtaining Correct Binding Free Energies in the Presence of Protein Conformational Change", vol. 3, no. 4, pp. 1231–1235, Jul. 2007, doi: 10.1021/ct700032n.
  27. J. Chen et al., "How Proximal Nucleobases Regulate the Catalytic Activity of G‑Quadruplex/Haemin DNAzymes", vol. 8, no. 12, pp. 11352–11361, Dec. 2018, doi: 10.1021/acscatal.8b03811.
  28. J. T. Hynes, "Chemical Reaction Dynamics in Solution", vol. 36, no. 1, pp. 573–597, Oct. 1985, doi: 10.1146/annurev.pc.36.100185.003041.