Drug Designing: Open Access
Open Access

Interaction, Dynamics and Stability Analysis of Some Minor Groove Binders with B-DNA Dodecamer 5' (CGCAAATTTGCG)-3'

Author info »

Introduction

DNA minor groove is the target for many anticancer and antitumor drugs. The forces that are involved in DNA-minor groove binding are electrostatic, van der Waal’s, hydrophobic, and hydrogen bonding [1]. Sequence specificity is important key tool for the drug-DNA interaction. Another important factor for the interaction is that the small molecule has a crescent shape which is complementary to the minor groove of DNA and thus minor groove binding is favored [2]. Frequently, minor groove shows selectivity towards AT region, several factors are responsible for this, first one is the electrostatic potential of AT-rich region is higher than that of GC-filled ones and another one is that the AT-rich grooves are narrower and deeper than GC ones [3]. DNA minor groove with alternating A and T, allows favorable van der Waal’s contacts between the drug and DNA instead of GC-rich region where the geometry of groove is altered by bulky amino groups of guanine bases [4,5].

In the current study, some different class of minor groove binder's (2,4-bis(4-amidinophenyl) furans and reversed diamidino 2,5-diarylfurans), labeled as mol-1, mol-2 mol-3 and mol-4 are taken, followed by a DNA sequence (PDB Id: 4AH0). Molecular modeling techniques such as geometry optimization, molecular docking and molecular dynamics simulation were used for the current research.

The ligands were first geometrically optimised in order for them to attain a minimum potential, as this would lead them to have electrostatically optimized chemical structure. Molecular Docking is simplest computational representation of DNA- ligand interactions [6]. After geometry optimization, ligands were docked to the selected DNA sequence for the identification of possible binding site followed by corresponding site binding energy was also obtained.

Various literatures have reported that, the minimum the site binding energy, the better the corresponding docked pose [7]. After geometry optimization and docking, the drug-DNA complexes were put to molecular dynamics simulation for 5ns and was observed that convergence is achieved within 5ns simulation. Obtained results not only predicted the conformational stability of the DNA with the ligands having had inserted between its base pairs but also gave a firm affirmation regarding the time dependence of the interaction and stability of drug-DNA complexes [8].

This article ends on a note that the three out of the four selected ligands had formed stable complexes with DNA; none of them had ruptured the DNA double helix over 5ns time scale simulation; interactions between drug and DNA was strong and significant; and helps in enriching the computational studies of drug-DNA interactions database and also paves ways for further sophisticated research viz, QM and QM/MM calculations [9,10]. The extension of such studies involving computational techniques is not only limited to drug design and molecular modelling but also to drug metabolism involving enzymes and their involved catalytic mechanisms [11].

Methods

System selection and preparation

The lead molecules (mol-1, mol-2, mol-3 and mol-4) were selected from literature [12]. Their chemical structures are shown in Figure 1. The DNA sequence; 4AH0 [13] was obtained from Protein Data Bank [14]. Its structural data including its specific nucleic acid bases sequences and PDB Id. is mentioned in Table 1 as obtained from Discovery Studio software package [15]. Water molecules and other residues, from selected DNA sequence were removed using UCSF Chimera before the start of docking calculations [16].

Figure 1: Chemical structures of the selected ligands.

Table 1: Nucleic Acid Report for 4AH0 as obtained from Discovery Studio Software.

Molecular geometry optimization

Using optimized ligand geometry before docking and other computational calculations is of utmost importance. The reason behind this is that geometry optimization would put the entire system in the lowest possible energy state; having least steric hindrances and charge. Further, the electrostatic charges also tend to get well distributed throughout the surface of the system leading to the electron transfer between definite pair of atoms; and eventually no distortions in the chemical structures are offered [17]. Here, the geometry optimization of selected ligands was carried out using Gaussian 09 software [18] using B3LYP hybrid functional at 6-31G** level of theory for their potentials to attain a local minima. These optimized ligands were then docked to the selected DNA sequence. Frequency calculation at the above level characterize the obtained point as minima on the potential energy surface in the optimize geometries.

Molecular docking studies

Molecular Docking is a widely used computational technique to predict the binding pocket of potent inhibiters in the vicinity of target molecules. The ability to determine the binding affinity of protein-ligand complexes is the key objective of computational docking [19]. Here, molecular docking was carried out using Autodock 4 software [20]. Docking calculations were set up using Lamarckian Genetic Algorithm (LGA). Gasteiger charges were added to each drug-DNA complex using Autodock Tools (ADT) before starting the docking simulations. A grid box, having various dimensions along the three coordinate axes, was prepared for each drug-DNA complex that enclosed the macromolecule. A 20 LGA run with a maximum cycle of 2500000 energy evaluations was performed for each of the drug-DNA complex. The docked pose with the least binding affinity was extracted and aligned with the receptor (DNA) for further analysis.

Molecular dynamics simulation

Molecular Dynamics is computational technique that is to study the dynamic characteristics, viz, protein folding, unwinding and other conformational changes including complex stability, over a specified period of time. Its applications have gained significant importance due to lack of experimental resources [21]. In the present work, GROMACS 5.0.4 software package [22] was used to carry out the molecular dynamics simulations. A total of four drug-DNA complexes were generated, after docking for molecular dynamics simulations viz, (4AH0+mol-1, 4AH0+mol-2, 4AH0+mol-3 and 4AH0+mol-4). These drug-DNA complexes were put to for 5000ps time scale simulation. There are many studies for the comparison of force fields for the nucleic acids but AMBER force fields [23] seems to be good for nucleic acid simulation due to the presence of specific topologies for the terminal nucleotides. Amber 03 force field already embedded in GROMACS software suite was used to generate the topology for the selected DNA sequence and antechamber module of AMBER program was used to generate topology of selected ligands through a python script: ‘acpype.py’ [24].The Ligand-DNA complex was solvated in a box of varying dimensions using TIP3P water model at 298K [25]. Sodium ions were then added to the solvated box containing the DNA-ligand complex by randomly replacing the water molecules in order to neutralize the system. Particle Mesh Ewald (PME) was used to handle long-range electrostatic interactions in periodic boundary conditions [26]. Energy minimization of the whole system was carried out in 25000 steps using Steepest Descent leap- Frog Integration Method followed by NVT ensemble equilibration at a constant temperature of 300K for 50s using Berendsen thermostat [27]. The system was then equilibrated with NPT ensemble at a constant pressure of 1atm in 25000 steps using steepest descent leap-frog integrator [27]. Particle Mesh Ewald (PME) was used to handle long-range electrostatic interactions in periodic boundary conditions and all the bonds involving hydrogen atoms were constrained using the LINCS algorithm [28]. Graphs were plotted using XMgrace software [29].

Natural bond orbital analysis

NBO analysis played the role of intermolecular orbital interaction in the complex, particularly charge transfer in the given wave function into localized form, corresponding to the one-center (lone pairs) and two-center (bonds) elements of the chemist's Lewis structure picture. In NBO analysis, the input atomic orbital basis set is transformed via natural atomic orbitals (NAOs) and natural hybrid orbitals (NHOs) into natural bond orbitals (NBOs). The NBOs obtained in this fashion correspond to the widely used Lewis picture, in which two-center bonds and lone pairs are localized. This is carried out by considering all possible interactions between filled donor NBOs and empty acceptor NBOs and estimating their energetic importance by second-order perturbation theory. For each donor NBO (i) and acceptor NBO (j), the hyper energy conjugated interaction (stabilization gained energy E⁽²⁾ in kcal/mol) associated with electron delocalization between donor NBOs and acceptor NBOs is estimated as

Where q_i is the orbital occupancy, ε_i, ε_j are diagonal elements and F_i,j is the off-diagonal NBO Fock matrix element [30].

Results and Discussion

Present study was aimed to identify new leads, targeting binding affinity, structural stability and applicability for antimicrobial ligands with DNA. The results obtained through various computational calculations are discussed and summarized as follows:

Molecular geometry optimization

The geometry of the lead molecules needs to be optimized before docking so that the repulsions between the nearest bonded atoms, between their angles and between their dihedrals could get minimized and the lead molecules would attain a state having least repulsion and hence owing to maximum stability and correspond to minimum energy [30]. The optimized geometries of the selected ligands are shown below in Figure 2.

Figure 2: Optimized structures of the selected ligands.

Molecular docking

The selected possess of ligands (mol-1, mol-2, mol-3 & mol-4) for antimicrobial tendencies were docked to 4AH0 in search of best docked posed complex leading to stability. The docking results, corresponding to the selected DNA sequences are summarized below in Table 2. Docking calculations revealed that 4AH0 with mol-4 had least binding energy of -12.39 kcal/mol. This indicates that since mol-4 had lesser binding energy, it owes to maximum stability of the formed complex. The docked pose corresponding to each drug-DNA complex is shown in Figure 3. These figures reveal that all the four ligands bonded themselves to the minor groove of each DNA sequence. Docking calculations also revealed that as bulky groups were attached to base structure of furan and then docked to DNA there were less no of hydrogen bonding interactions; owing to more steric hindrances and consequently the fact that number of hydrogen bonding interactions decrease on attaching bulky groups complement the obtained results. These hydrogen bonding interactions were obtained from Discovery Studio Visualizer [31].

Figure 3: Best docked posed complexes for 4AH0 as obtained from UCSF Chimera.

Table 2: Various docking results obtained for 4AH0DNA sequence.

Following Figures 4 and 5 represent the atomistic binding sites for each of the drug-DNA complexes and corresponding H-donor/ acceptor clouded regions at binding sites, due care while taking images regarding the same has been taken in order to have the same pose for both, binding site and H-donor/acceptor cloud, respectively. Further Table 3 represents the donor and the acceptor residues involved in the formation of hydrogen bond between the DNA and ligands.

Figure 4: H-bonds in best docked posed complexes for 4AH0 as obtained from Discovery Studio Visualizer.

Figure 5: H-bond donor and acceptor regions in best docked posed complexes for 4AH0 as obtained from Discovery Studio Visualizer.

Table 3: The donor and the acceptor amino acid residues involved in the formation of hydrogen bond between the DNA and ligand atoms.

Molecular dynamics simulation

Structural stability of biomolecules under dynamical conditions over a specifically mimicked environment for a pre-defined period of time can be studied via molecular dynamics simulations. Such studies hold significant importance in the structure and dynamics of biomolecules owing to less experimental costs. Various parameters were studied and analyzed for their contribution to stability analysis of the drug-DNA complexes, and are discussed as follows:

Variation in radius of gyration: In order to understand the dynamic stability and compactness of DNA-ligand complexes, Radii of Gyration values are determined. Clearly, the avg. radius of gyration lies between 1.25 nm to 1.50 nm resp. The collective Radius of Gyration variations for each drug-DNA pairs are shown in Figure 6. The Radius of Gyration data reveals that mol-2, mol-3 and mol-4 with 4AH0 sequence form themost stable complexes whereas for mol-1, the radius remains intact till ~3500ps but after that time limit there are frequent perturbations in the conformations of the complex owing to minimized stability.

Figure 6: Variations for Radius of Gyration.

Root mean square deviation: The RMSD is treated as the measure of the conformational stability of biological macromolecules. The plots for RMSD of all the drug- DNA complexes are represented in Figure 7. From the variations shown below, the RMSD for corresponding to mol-2, mol-3 and mol-4 with 4AH0 almost lie between 0 to 0.5 nm representing least or no deviations. However,4AH0+mol-1, shows deviations around 3400ps, 4000ps and 4500ps and hence confirms for having had formed least stable complex with considerable deviations in the DNA double helix.

Figure 7: RMSD for Drug-DNA Complexes.

Root mean square fluctuation: RMSF records the fluctuation of each amino acid base pair including the fluctuations in the flexible regions within the nucleic acid during the course of the MD simulation.Well-structured regions and loosely bound regions in DNA strands are distinguished by low and high root mean square fluctuations values respectively. For small proteins, a fluctuation between 1~3 Å is acceptable. The graphs shown below in Figure 8 suggest that, 4AH0+mol-2, 4AH0+mol3 and 4AH0+mol-4 had the least fluctuations, and especially around and beyond residue number 800 among and hence owe towards the stability of the complexes.Here this range lies between 0~0.5Å. Whereas,4AH0+mol-1 shows maximum deviations as seen from the corresponding graph.

Figure 8: RMSF for Drug-DNA Complexes.

Natural bond orbital analysis: Natural Bond Orbital (NBO) analysis played a crucial role for the calculation of intra and intermolecular orbital bonding and interaction in the optimized electronic chemical structure, particularly charge transfer and stabilization gain energy. Geometry, of the NBO analysis were optimized at the level of B3LYP/6-31G (d) method using Gaussian 09 program package (Figure 9).

Figure 9: Optimized geometries are showing the molecular orbital with maximum electron density.

The NBO analysis is a suitable source of the charge transfer donner to receptor in the orbital-based. The total hyper energy conjugated interaction (stabilization gained energy) from NBO donor to NBO acceptors were calculated by the second-order perturbative theory. The estimates of donor-acceptor interactions in the NBO basis and presented in the Table supporting information (SI) Table S1-S4. This analysis clearly indicates that the stabilization gained energy is a few times higher in the total stabilization gain energy [32].

Conclusion

Modeling techniques when applied in the studies of drugs and associated targets have always proven them to be key players in pharmacy industries. These techniques not only not only complement experiments but also help in the improvement of existing computational tools.

The current research article addressed a problem regarding interaction and stability of drugs in the vicinity of DNA, and our findings successfully resolve this claim. The optimized drugs structures when docked to selected DNA generated best docked pose by finding the most stable binding site, these docked complexes when subjected to molecular dynamics for time scale evaluation of their stability yielded similar results that complemented the docking results and previously reported datas. Our evaluated parameters viz, hydrogen bonding, Radius of Gyration, RMSD and RMSF were sufficient to support our claim. We successfully traced out the most stable complexes, 4AH0+mol-2, 4AH0+mol-3 and 4AH0+mol-4. This clearly raises numerous questions to be well addressed regarding their human bioassays, but we can say that in vivo experiments are solely responsible for any sure and certain conclusions. However, computational techniques have a very bright future in this field. This not only keeps the doors open for further research but also leaves us with strong and affirmative methodologies. The NBO analysis is provided the detailed insight into the type of hybridization and the nature of bonding in the particularly charge transfer and stabilization gain energy.

Conflict of Interest

The authors declare no conflicts in their mutual interests, in any manner whatsoever.

Acknowledgements

Anwesh Pandey acknowledges to the University Grants Commission, India for financial assistance during the course of this research work and Suresh Kumar is also grateful to the University Grants Commission for financial support through the Dr. D. S. Kothari Postdoctoral Fellowship (No.F.4-2/2006 (BSR)/PH/18-19/0022).

References

1. Watson JD, Crick FH. Molecular structure of nucleic acids: a structure for deoxyribose nucleic acid. Nature. 1953; 171(4356):737-738.
2. Pindur U, Jansen M, Lemster T. Advances in DNA-ligands with groove binding, intercalating and/or alkylating activity: chemistry, DNA-binding and biology. Curr Med Chem. 2005; 12(24):2805-2847.
3. Yadava U, Yadav SK, Yadav RK. Investigations on bisamidine derivatives as novel minor groove binders with the dodecamer 5′(CGCGAATTCGCG) 3′. J Mol Liq. 2019; 280:135-152.
4. Wilson WD, Tanious FA, Mathis A, Tevis D, Hall JE, Boykin DW. Antiparasitic compounds that target DNA. Biochimie. 2008; 90(7):999-1014.
5. Yadava U, Yadav SK, Yadav RK. Electronic structure, vibrational assignments and simulation studies with A/T rich DNA duplex of an aromatic bis-amidine derivative. DNA repair. 2017; 60:9-17.
6. Pandey A, Mishra R, Shukla A, Yadav AK, Kumar D. Proceeding of international symposium on advances in functional and biological materials (ISAFBM-2019). ISAFBM Conference Proceedings 2019.
7. Srivastava HK, Chourasia M, Kumar D, Sastry GN. Comparison of computational methods to model DNA minor groove binders. J Chem Inf Model. 2011; 51(3):558-571.
8. Pandey A, Yadav R, Shukla A, Yadav AK. Unveiling the Antimicrobial Activities of Dicationic Carbazoles and Related Analogs Through Computational Docking. Adv Sci Eng Med. 2020; 12(1):40-44.
9. Senn HM, Thiel W. QM/MM methods for biomolecular systems. Angew Chem Int Ed. 2009; 48(7):1198-1229.
10. Shaik S, Kumar D, de Visser SP, Altun A, Thiel W. Theoretical perspective on the structure and mechanism of cytochrome P450 enzymes. Chem Rev. 2005; 105(6):2279-2328.
11. Yadav R, Pandey A, Awasthi N, Shukla A. Molecular Docking Studies of Enzyme Binding Drugs on Family of Cytochrome P450. Adv Sci Eng Med. 2020; 12(1):83-87.
12. Tidwell RR, Boykin DW. Dicationic DNA minor groove binders as antimicrobial agents. Small molecule DNA and RNA binders: from synthesis to nucleic acid complexes. J Braz Chem Soc. 2002; 27:414-460.
13. Balendiran K, Rao ST, Sekharudu CY, Zon G, Sundaralingam M. X-ray structures of the B-DNA dodecamer d (CGCGTTAACGCG) with an inverted central tetranucleotide and its netropsin complex. Acta Crystallogr D Biol Crystallogr. 1995; 51(2):190-198.
14. Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, et al. The protein data bank. Nucleic Acids Res. 2000; 28(1):235-242.
15. Pettersen EF, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC, et al. UCSF Chimera-a visualization system for exploratory research and analysis. J Comput Chem. 2004; 25(13):1605-1612.
16. Dubey DK, Chaubey KA, Parveen A, Ojha PR. Comparative study of inhibition of drug potencies of c-Abl human kinase inhibitors: a computational and molecular docking study. J Biophys Struct Biol. 2010; 2(4):47-54.
17. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, et al. Electronic supplementary material (ESI) for chemical science. R Soc Chem. 2016.
18. Mishra R, Kumar A, Chandra R, Kumar D. A review on theoretical studies of various types of Drug-DNA Interaction. Int J Sci Technol Soc. 2017; 3:11-27.
19. Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS, et al. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J Comput Chem. 2009; 30(16):2785-2791.
20. Zhao H, Caflisch A. Molecular dynamics in drug design. Eur J Med Chem. 2015; 91:4-14.
21. Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, Lindahl E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX. 2015; 1:19-25.
22. Ponder JW, Case DA. Force fields for protein simulations. Adv Prot Chem. 2003; 66:27-85.
23. Da Silva AW, Vranken WF. ACPYPE-Antechamber python parser interface. BMC Res Notes 2012; 5(1):367.
24. Mark P, Nilsson L. Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K. J Phys Chem A. 2001; 105(43):9954-9960.
25. Darden T, York D, Pedersen L. Particle mesh Ewald: An N⋅ log (N) method for Ewald sums in large systems. J Chem Phys. 1993; 98(12):10089-10092.
26. Ryckaert JP, Ciccotti G, Berendsen HJ. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J Comput Phys. 1977; 23(3):327-341.
27. Hess B, Bekker H, Berendsen HJ, Fraaije JG. LINCS: a linear constraint solver for molecular simulations. J Comput Chem. 1997; 18(12):1463-1472.
28. Turner PJ. XMGRACE, Version 5.1.19. Center for coastal and land-margin research, oregon graduate institute of science and technology, Beaverton, OR. 2005.
29. Gangadharan RP, Krishnan SS. Natural bond orbital (nbo) population analysis of 1-azanapthalene-8-ol. Acta Phys Pol A. 2014; 125(1).
30. Shastri R, Awasthi N, Kumar D, Yadav AK, Roy D, Goutam SP, et al. A density functional theory study on structural stability and electronic properties of co x o y (x+ y= 4-12) nanoclusters. Adv Sci Eng Med. 2018; 10(12):814-818.
31. Rather AM, Shome A, Kumar S, Bhunia BK, Mandal BB, Srivastava HK, et al. Alkali metal-ion assisted Michael addition reaction in controlled tailoring of topography in a superhydrophobic polymeric monolith. J Mater Chem A. 2018; 6(35):17019-17031.