The molecular mechanism of ligand unbinding from the human telomeric G-quadruplex by steered molecular dynamics and umbrella sampling simulations

The interactions of Pu22 G-quadruplex, derived from c-MYC promoter sequence, with antitumor acridine derivatives-An NMR/MD combined study

Unsymmetrical bisacridines (UAs) represent a novel class of anticancer agents that exhibit significant antitumor activity against a wide range of cancer cell lines and solid tumors in vivo. UAs consist of two different acridine-based ring systems, which are connected by an aminoalkyl linker. Recent studies have demonstrated that UAs can suppress the c-MYC protooncogene, which is overexpressed in many tumor types. As a proposed molecular basis for this activity, UAs have been suggested to stabilize the G-quadruplex structure formed within the promoter region of c-MYC. In this study, we performed spectroscopic and computational analyses to investigate the stereochemistry of the c-MYC NHE III1 representative G-quadruplex, codenamed Pu22, in complex with two promising bisacridines, C-2045 and C-2053, as well as their monomeric counterparts, C-1311 and C-1748. C-1311 formed a well-defined 1:2 mol/mol DNA:ligand non-covalent adduct, whose solution structure was determined via 2D NMR. In contrast, C-1748 displayed weak and nonspecific interactions with the Pu22 G-quadruplex. Finally, the Pu22:UA complexes were examined using a combination of NMR and molecular modeling approaches, including umbrella sampling simulations. These results provide insights into the interaction mechanisms of UAs with G-quadruplex structures and highlight their potential as therapeutic agents targeting c-MYC.

Post-Docking Refinement of Peptide or Protein-RNA Complexes Using Thermal Titration Molecular Dynamics (TTMD): A Stability Insight

RNA-protein interactions drive and regulate fundamental cellular processes like transcription and translation. Despite being still limited, the growing body of structural data significantly contributes to the characterization of these interactions. However, RNA complexes involving proteins or peptides are not always available due to the structural determination challenges that this biopolymer entails. Consequently, modeling approaches like molecular docking are exploited to generate complexes relevant to structural and pharmaceutical purposes, including analysis of putative drug targets. Docking methods, despite their widespread adoption, are often hindered by limitations in scoring accuracy, which affects the ranking of the generated poses. Postdocking refining methods, including molecular dynamics (MD) approaches, have been developed to tackle this issue. Thermal Titration Molecular Dynamics (TTMD) is an enhanced sampling molecular dynamics technique that has been previously effectively applied to refine protein or RNA-small-molecule docking poses. This study presents the first application of TTMD to RNA-peptide complexes, validating this method on more complex systems and extending its applicability domain. Our findings showcase the capability of this technique to refine peptide-RNA docking poses, correctly identifying native binding modes among decoys for different pharmaceutically relevant targets.

Microscopic Understanding of the Conformational Stability of the Aggregated Nonamyloid β Components of α-Synuclein

Self-association of α-synuclein peptides into oligomeric species and ordered amyloid fibrils is associated with Parkinson's disease, a progressive neurodegenerative disorder. In particular, the peptide domain formed between the residues Glu-61 (or E61) and Val-95 (or V95) of α-synuclein, typically termed the "nonamyloid β component" (NAC), is known to play critical roles in forming aggregated structures. In this work, we have employed molecular dynamics simulations to explore the conformational properties and relative stabilities of aggregated protofilaments of different orders, namely, tetramer (P(4)), hexamer (P(6)), octamer (P(8)), decamer (P(10)), dodecamer (P(12)), and tetradecamer (P(14)), formed by the NAC domains of α-synuclein. Besides, center-of-mass pulling and umbrella sampling simulation methods have also been employed to characterize the mechanistic pathway of peptide association/dissociation and the corresponding free energy profiles. Structural analysis showed that the disordered C-terminal loop and the central core regions of the peptide units lead to more flexible and distorted structures of the lower order protofilaments (P(4) and P(6)) as compared to the higher order ones. Interestingly, our calculation shows the presence of multiple distinctly populated conformational states for the lower order protofilament P(4), which may drive the oligomerization process along multiple pathways to form different polymorphic α-synuclein fibrillar structures. It is further observed that the nonpolar interaction between the peptides and the corresponding nonpolar solvation free energy play a dominant role in stabilizing the aggregated protofilaments. Importantly, our result showed that reduced cooperativity during the binding of a peptide unit beyond a critical size of the protofilament (P(12)) leads to less favorable binding free energy of a peptide.

Deciphering the competitive inhibition of dihydropteroate synthase by 8 marcaptoguanine analogs: enhanced potency in phenylsulfonyl fragments

The emergence of sulfa-drug resistance and reduced efficacy of pterin-based analogs towards Dihydropteroate synthase (DHPS) inhibition dictate a pressing need of developing novel antimicrobial agents for immune-compromised patients. Recently, a series of 8-Marcaptoguanin (8-MG) derivatives synthesized for 6-Hydroxymethyl-7,8-dihydropterin pyrophosphokinase (experimental K_D ∼ 100-.0.36) showed remarkable homology with the pteroic-acid and serve as a template for product antagonism in DHPS. The present work integrates ligand-based drug discovery techniques with structure-based docking, enhanced MD simulation, and MM/PBSA techniques to demonstrate the essential features of 8-MG analogs which make it a potent inhibitor for DHPS. The delicate balance in hydrophilic, hydrophobic substitutions on the 8-MG core is the crucial signature for DHPS inhibition. It is found that the dynamic interactions of active compounds are mainly dominated by consistent hydrogen bonding network with Asp 96, Asn 115, Asp 185, Ser 222, Arg 255 and π-π stacking, π-cation interactions with Phe 190, Lys 221. Further, two new 8-MG compounds containing N-phenylacetamide (compound S1, ΔG_bind-eff = -62.03 kJ/mol) and phenylsulfonyl (compound S3, ΔG_bind-eff = -71.29 kJ/mol) fragments were found to be the most potent inhibitor of DHPS, which stabilize the flexible pABA binding loop, thereby increasing their binding affinity. MM/PBSA calculation shows electrostatic energy contribution to be the principal component in stabilizing the inhibitors in the binding pocket. This fact is further confirmed by the higher energy barrier obtained in umbrella sampling for this class of inhibitors.

The Molecular Mechanism of Human Voltage-Dependent Anion Channel 1 Blockade by the Metallofullerenol Gd@C82(OH)22: An In Silico Study

The endohedral metallofullerenol Gd@C₈₂(OH)₂₂ has been identified as a possible antineoplastic agent that can inhibit both the growth and metastasis of cancer cells. Despite these potentially important effects, our understanding of the interactions between Gd@C₈₂(OH)₂₂ and biomacromolecules remains incomplete. Here, we study the interaction between Gd@C₈₂(OH)₂₂ and the human voltage-dependent anion channel 1 (hVDAC1), the most abundant porin embedded in the mitochondrial outer membrane (MOM), and a potential druggable target for novel anticancer therapeutics. Using in silico approaches, we observe that Gd@C₈₂(OH)₂₂ molecules can permeate and form stable interactions with the pore of hVDAC1. Further, this penetration can occur from either side of the MOM to elicit blockage of the pore. The binding between Gd@C₈₂(OH)₂₂ and hVDAC1 is largely driven by long-range electrostatic interactions. Analysis of the binding free energies indicates that it is thermodynamically more favorable for Gd@C₈₂(OH)₂₂ to bind to the hVDAC1 pore when it enters the channel from inside the membrane rather than from the cytoplasmic side of the protein. Multiple factors contribute to the preferential penetration, including the surface electrostatic landscape of hVDAC1 and the unique physicochemical properties of Gd@C₈₂(OH)₂₂. Our findings provide insights into the potential molecular interactions of macromolecular biological systems with the Gd@C₈₂(OH)₂₂ nanodrug.

Efficient Estimation of Absolute Binding Free Energy for a Homeodomain-DNA Complex from Nonequilibrium Pulling Simulations

Estimation of binding free energies is one of the central aims of simulations of biomolecular complexes. We explore the accuracy and efficiency of setups based on nonequilibrium pulling simulations applied to the estimation of binding affinities of DNA-binding proteins. Absolute binding free energies are calculated over a range of temperatures and compared to results obtained previously using an equilibrium method. We show that realistic binding affinities can be obtained with the presented nonequilibrium approach, which also entails lower computational requirements. Errors of the binding free energy estimates are investigated and are shown to be comparable to those observed previously. Bounds are provided on the convergence of the errors with respect to the number of pulling simulations performed and with respect to the applied pull rate.

Enhanced sampling and free energy calculations for protein simulations

Molecular dynamics simulation is a powerful computational technique to study biomolecular systems, which complements experiments by providing insights into the structural dynamics relevant to biological functions at atomic scale. It can also be used to calculate the free energy landscapes of the conformational transitions to better understand the functions of the biomolecules. However, the sampling of biomolecular configurations is limited by the free energy barriers that need to be overcome, leading to considerable gaps between the timescales reached by MD simulation and those governing biological processes. To address this issue, many enhanced sampling methodologies have been developed to increase the sampling efficiency of molecular dynamics simulations and free energy calculations. Usually, enhanced sampling algorithms can be classified into methods based on collective variables (CV-based) and approaches which do not require predefined CVs (CV-free). In this chapter, the theoretical basis of free energy estimation is briefly reviewed first, followed by the reviews of the most common CV-based and CV-free methods including the presentation of some examples and recent developments. Finally, the combination of different enhanced sampling methods is discussed.

Why do G-quadruplexes dimerize through the 5'-ends? Driving forces for G4 DNA dimerization examined in atomic detail

G-quadruplexes (G4) are secondary structures formed by guanine-rich nucleic acid sequences and shown to exist in living cells where they participate in regulation of gene expression and chromosome maintenance. G-quadruplexes with solvent-exposed guanine tetrads show the tendency to associate together through cofacial stacking, which may be important for packaging of G4-forming sequences and allows for the design of higher-order G4 DNA structures. To understand the molecular driving forces for G4 association, here, we study the binding interaction between two parallel-stranded G-quadruplexes using all-atom molecular dynamics simulations. The predicted dimerization free energies show that direct binding through the 5’-G-tetrads is the most preferred of all possible end-to-end stacking orientations, consistently with all available experimental data. Decomposition of dimerization enthalpies in combination with simulations at varying ionic strength further indicate that the observed orientational preferences arise from a fine balance between the electrostatic repulsion of the sugar-phosphate backbones and favorable counterion binding at the dimeric interface. We also demonstrate how these molecular-scale findings can be used to devise means of controlling G4 dimerization equilibrium, e.g., by altering salt concentration and using G4-targeted ligands.

Native DNA usually folds to form the canonical double helix, however, under certain conditions, it can also fold into other secondary structures. Some of the most interesting ones are G-quadruplexes (G4)—compact DNA structures in which guanines assemble into multilayered tetrads, and whose formation has been reported at the ends of linear chromosomes (telomeres) and at different regulatory regions of the genome. Although structural and basic energetic properties, as well as some biological functions of G-quadruplexes are quite well understood, not much is known about their propensity to form agregated structures. A very high density of G-quadruplexes at telomeres along with their large exposed planar surfaces indeed favor G4 aggregation through end-to-end stacking, which might be important for the protection of telomeres and DNA packaging. In this research, using computer simulations, we provide insight into molecular origins of stability of the higher-order G-quadruplexes and explain in structural and energetic terms a strong preference for one particular end-to-end stacking orientation. Based on the recognized aggregation driving forces, we also suggest methods for controling the aggregation preferences openining up new opportunities for designing oligomeric G-quadruplexes.

Surface Topography Effects of Globular Biomolecules on Hydration Water

Hydration water serves as a microscopic manifestation of structural stability and functions of biomolecules. To develop bio-nanomaterials in applications, it is important to study how the surface topography and heterogeneity of biomolecules result in their diversity of the hydration dynamics and energetics. We here performed molecular dynamics simulations combined with the steered molecular dynamics and umbrella sampling to investigate the dynamics and escape process associated with the free energy change of water molecules close to a globular biomolecule, i.e., hemoglobin (Hb) and G-quadruplex DNA (GDNA). The residence time, power of long-time tail, and dipole relaxation time were found to display drastic changes within the averaged hydration shell of 3.0-5.0 Å. Compared with bulk water, in the inner hydration shell, the water dipole moment displays a slower relaxation process and is more oriented toward GDNA than toward Hb, forming a hedgehog-like structure when it surrounds GDNA. In particular, a spine water structure is observed in the GDNA narrow groove. The water isotope effect not only prolongs the dynamic time scales of libration motion in the inner hydration shell and the dipole relaxation processes in the bulk but also strengthens the DNA spine water structure. The potential of the mean force profile reflects the integrity of the hydration shell structure and enables us to obtain detailed insights into the structures formed by water, such as the caged H-bond network and the edge bridge structures; it also reveals that local hydration shell free energy (LHSFE) depends on H-bond rupture processes and ranges from 0.2 to 4.2 kcal/mol. Our results demonstrate that the surface topography of a biomolecule influences the integrity of the hydration shell structure and LHSFE. Our studies are able to identify various further applications in the areas of microfluid devices and nano-dewetting on bioinspired surfaces.

Potential Mean Force from Umbrella Sampling Simulations: What Can We Learn and What Is Missed?

Changes in free energy provide valuable information for molecular recognition, including both ligand-receptor binding thermodynamics and kinetics. Umbrella sampling (US), a widely used free energy calculation method, has long been used to explore the dissociation process of ligand-receptor systems and compute binding free energy. In existing publications, the binding free energy computed from the potential of mean force (PMF) with US simulation mostly yielded "ball park" values with experimental data. However, the computed PMF values are highly influenced by factors such as initial conformations and/or trajectories provided, the reaction coordinate, and sampling of conformational space in each US window. These critical factors have rarely been carefully studied. Here we used US to study the guest aspirin and 1-butanol dissociation processes of β-cyclodextrin (β-CD) and an inhibitor SB2 dissociation from a p38α mitogen-activated protein kinase (MAPK) complex. For β-CD, we used three different β-CD conformations to generate the dissociation path with US windows. For p38α, we generated the dissociation pathway by using accelerated molecular dynamics followed by conformational relaxing with short conventional MD, steered MD, and manual pulling. We found that, even for small β-CD complexes, different β-CD conformations altered the height of the PMF, but the pattern of PMF was not affected if the MD sampling in each US window was well-converged. Because changing the macrocyclic ring conformation needs to rotate dihedral angles in the ring, a bound ligand largely restrains the motion of cyclodextrin. Therefore, once a guest is in the binding site, cyclodextrin cannot freely change its initial conformation, resulting in different absolute heights of the PMF, which cannot be overcome by running excessively long MD simulations for each US window. Moreover, if the US simulations were not converged, the important barrier and minimum were missed. For ligand-protein systems, our studies also suggest that the dissociation trajectories modeled by an enhanced sampling method must maintain a natural molecular movement to avoid biased PMF plots when using US simulations.

Tuning the ring strain effect in acridine derivatives on binding affinity with G-quadruplex-DNA: A computational and experimental study

Search for inhibitors to stabilize the telomeric G-quadruplex in order to deter telomerase activity is an active area of research. Inhibitors play an important role to initiate the tumor cell mortalization process. This work reports for the first time of acridine derivative with four membered ammonium rings at the side chain to surpass the binding ability against BRACO-19 with G-quadruplex-DNA. It is known in the literature that acridine based molecule BRACO-19 can effectively bind with G-quadruplex-DNA. The computational study performed in this study revealed that the binding ability of acridine based molecule can be augmented with subtle variation in the molecular structure of the drug like candidates. Steered molecular dynamics (SMD) performed with the acridine derivatives and G-quadruplex DNA showed the importance of ring strain to the side chain of those ligand molecules. The rupture force analysis, hydrogen bonding interactions and the calculated free energies in MM-PBSA method suggest that ligand 3 is superior than that of BRACO-19. The synthesized ligand 3 and BRACO-19 showed the binding constants obtained from ITC measurements are 4 × 10⁶ mol^-1 and 2.6 × 10⁶, which corroborates the computational findings.

Binding modes and pathway of RHPS4 to human telomeric G-quadruplex and duplex DNA probed by all-atom molecular dynamics simulations with explicit solvent

RHPS4, a potent binder to human telomeric DNA G-quadruplex, shows high efficacy in tumor cell growth inhibition. However, it's preferential binding to DNA G-quadruplex over DNA duplex (about 10 fold) remains to be improved toward its clinical application. A high resolution structure of the single-stranded telomeric DNA G-quadruplexes, or B-DNA duplex, in complex with RHPS4 is not available yet, and the binding nature of this ligand to these DNA forms remains to be elusive. In this study, we carried out 40 μs molecular dynamics binding simulations with a free ligand to decipher the binding pathway of RHPS4 to a DNA duplex and three G-quadruplex folders (parallel, antiparallel and hybrid) of the human telomeric DNA sequence. The most stable binding mode identified for the duplex, parallel, antiparallel and hybrid G-quadruplexes is an intercalation, bottom stacking, top intercalation and bottom intercalation mode, respectively. The intercalation mode with similar binding strength to both the duplex and the G-quadruplexes, explains the lack of binding selectivity of RHPS4 to the G-quadruplex form. Therefore, a ligand modification that destabilizes the duplex intercalation mode but stabilizes the G-quadruplex intercalation mode will improve the binding selectivity toward G-quadruplex. The intercalation mode of RHPS4 to both the duplex and the antiparallel and the hybrid G-quadruplex follows a base flipping-insertion mechanism rather than an open-insertion mechanism. The groove binding, the side binding and the intercalation with flipping out of base were observed to be intermediate states before the full intercalation state with paired bases.

Resolving the Ligand-Binding Specificity in c-MYC G-Quadruplex DNA: Absolute Binding Free Energy Calculations and SPR Experiment

We report the absolute binding free energy calculation and surface plasmon resonance (SPR) experiment for ligand binding with the c-MYC G-quadruplex DNA. The unimolecular parallel DNA G-quadruplex formed in nuclease hypersensitivity element III₁ of the c-MYC gene promoter regulates the c-MYC transcription and is recognized as an emerging drug target for cancer therapy. Quindoline derivatives have been shown to stabilize the G-quadruplex and inhibit the c-MYC expression in cancer cells. NMR revealed two binding sites located at the 5' and 3' termini of the G-quadruplex. Questions about which site is more favored and the basis for the ligand-induced binding site formation remain unresolved. Here, we employ two absolute binding free energy methods, the double decoupling and the potential of mean force methods, to dissect the ligand-binding specificity in the c-MYC G-quadruplex. The calculated absolute binding free energies are in general agreement with the SPR result and suggest that quindoline has a slight preference for the 5' site. The flanking residues around the two sites undergo significant reorganization as the ligand unbinds, which provides evidence for ligand-induced binding pocket formation. The results help interpret experimental data and inform rational design of small molecules targeting the c-MYC G-quadruplex.

Probing the Binding Pathway of BRACO19 to a Parallel-Stranded Human Telomeric G-Quadruplex Using Molecular Dynamics Binding Simulation with AMBER DNA OL15 and Ligand GAFF2 Force Fields

Human telomeric DNA G-quadruplex has been identified as a good therapeutic target in cancer treatment. G-quadruplex-specific ligands that stabilize the G-quadruplex have great potential to be developed as anticancer agents. Two crystal structures (an apo form of parallel stranded human telomeric G-quadruplex and its holo form in complex with BRACO19, a potent G-quadruplex ligand) have been solved, yet the binding mechanism and pathway remain elusive. In this study, we simulated the binding of a free BRACO19 molecule to the apo form of the G-quadruplex using the latest AMBER DNA (OL15) and ligand (GAFF2) force fields. Three binding modes have been identified: top stacking, bottom intercalation, and groove binding. Bottom intercalation (51% of the population) resembles the bottom binding pose in the complex crystal structure very well. The groove binding mode is less stable than the bottom binding mode and is likely to be an intermediate state leading to the bottom binding mode. A flip-insertion mechanism was observed in the bottom intercalation mode, during which flipping of the bases outward makes space for ligand insertion, after which the bases flip back to increase the stability of the complex. In addition to reproducing the base-flipping behavior for some loop residues upon ligand binding, the direct alignment type of the ATAT-tetrad was observed in our simulations for the first time. These successes provide initial support for using this combination of the OL15 and GAFF2 force fields to study quadruplex-ligand interactions.

High-resolution three-dimensional NMR structure of the KRAS proto-oncogene promoter reveals key features of a G-quadruplex involved in transcriptional regulation

Non-canonical base pairing within guanine-rich DNA and RNA sequences can produce G-quartets, whose stacking leads to the formation of a G-quadruplex (G4). G4s can coexist with canonical duplex DNA in the human genome and have been suggested to suppress gene transcription, and much attention has therefore focused on studying G4s in promotor regions of disease-related genes. For example, the human KRAS proto-oncogene contains a nuclease-hypersensitive element located upstream of the major transcription start site. The KRAS nuclease-hypersensitive element (NHE) region contains a G-rich element (22RT; 5'-AGGGCGGTGTGGGAATAGGGAA-3') and encompasses a Myc-associated zinc finger-binding site that regulates KRAS transcription. The NEH region therefore has been proposed as a target for new drugs that control KRAS transcription, which requires detailed knowledge of the NHE structure. In this study, we report a high-resolution NMR structure of the G-rich element within the KRAS NHE. We found that the G-rich element forms a parallel structure with three G-quartets connected by a four-nucleotide loop and two short one-nucleotide double-chain reversal loops. In addition, a thymine bulge is found between G8 and G9. The loops of different lengths and the presence of a bulge between the G-quartets are structural elements that potentially can be targeted by small chemical ligands that would further stabilize the structure and interfere or block transcriptional regulators such as Myc-associated zinc finger from accessing their binding sites on the KRAS promoter. In conclusion, our work suggests a possible new route for the development of anticancer agents that could suppress KRAS expression.

Product Release Pathways in Human and Plasmodium falciparum Phosphoribosyltransferase

Atomistic molecular dynamics (MD) simulations coupled with the metadynamics technique were carried out to delineate the product (PPi.2Mg and IMP) release mechanisms from the active site of both human (Hs) and Plasmodium falciparum (Pf) hypoxanthine-guanine-(xanthine) phosphoribosyltransferase (HG(X)PRT). An early movement of PPi.2Mg from its binding site has been observed. The swinging motion of the Asp side chain (D134/D145) in the binding pocket facilitates the detachment of IMP, which triggers the opening of flexible loop II, the gateway to the bulk solvent. In PfHGXPRT, PPi.2Mg and IMP are seen to be released via the same path in all of the biased MD simulations. In HsHGPRT too, the product molecules follow similar routes from the active site; however, an alternate but minor escape route for PPi.2Mg has been observed in the human enzyme. Tyr 104 and Phe 186 in HsHGPRT and Tyr 116 and Phe 197 in PfHGXPRT are the key residues that mediate the release of IMP, whereas the motion of PPi.2Mg away from the reaction center is guided by the negatively charged Asp and Glu and a few positively charged residues (Lys and Arg) that line the product release channels. Mutations of a few key residues present in loop II of Trypanosoma cruzi (Tc) HGPRT have been shown to reduce the catalytic efficiency of the enzyme. Herein, in silico mutation of corresponding residues in loop II of HsHGPRT and PfHGXPRT resulted in partial opening of the flexible loop (loop II), thus exposing the active site to bulk water, which offers a rationale for the reduced catalytic activity of these two mutant enzymes. Investigations of the product release from these HsHGPRT and PfHGXPRT mutants delineate the role of these important residues in the enzymatic turnover.