Novel G-quadruplex stabilizing agents: in-silico approach and dynamics

The stabilization of overhang G-rich repetitive DNA units at the 3′-end of telomeres, which are well known to form functionally important G-quadruplex structures, is a current goal in designing novel anticancer drugs. In the present study, we have undertaken an in silico approach by molecular docking using a small molecule library to find potential G-quadruplex stabilizing agents. Two molecules, A, [N′1-imino(2-pyridyl)methyl-3,4,5-trimethoxybenzene-1-carbohydrazide] and B, [(3-[4-({[3-({4-[(2cyanoethyl)(methyl)amino]benzylidene}amino)propyl]imino}methyl)(methyl) anilino]propanenitrile)], that had good docking scores have been investigated for interaction with G-quadruplexes in a Molecular Dynamics simulation study. Fluorescence spectroscopy of G-quadruplexes bound to the screened molecules A and B was used to experimentally validate the theoretical results. The binding of ligands A and B to G-quadruplexes resulted in blue shifts of 10–18 nm, respectively, in the fluorescence emission spectra of the G-quadruplexes, demonstrating that both molecules bind to the G-face of the quadruplex. The same experiment was performed for the complexation of these small molecules with a G-rich DNA duplex, . Interestingly, no blue shift was observed in the fluorescence emission spectra of the DNA duplex in the presence of these small molecules. Thus, these findings indicated that these ligands very selectively bind to G-quadruplexes instead of the duplex DNA. In addition, a one-dimensional water ligand observed via a gradient spectroscopy Nuclear Magnetic Resonance (NMR) experiment showed that both molecules bound to the 23-mer G-quadruplex DNA. The molecular properties of the ligand–quadruplex complex have been analyzed with the help of the Adaptive Poisson-Boltzmann Solver, revealing that electrostatics govern the binding of the small molecules to G-quadruplexes. Both molecules were investigated in detail using solvation free energy calculations and Absorption, Distribution, Metabolism, Elimination and Toxicity (ADMET) predictions, which provide insight into lead optimization for designing G-quadruplex stabilizing agents; therefore, these molecules have potential as new therapeutic agents.


Introduction
The 3′-end region of chromosomes consists of almost 200 nucleotides composed of tandem repeats of hexanucleotide (TTAGGG) n and ends with an overhanging single strand. Telomeres are the main location in which these crucial functional secondary structures of DNA form (Cech, 2004). Telomeres have a vital role in life because of the various functions of telomere-associated proteins, their attachment to the nuclear matrix, and their higher order chromatin structures (Campbell, 2012;Collado, Blasco, & Serrano, 2007). Conservation of telomeric length is an important biological condition for cell growth (Engelhardt & Finke, 2001;Lange, 2009). Telomerase is the enzyme that functionalizes the addition of hexanucleotide repeats of TTAGGG to the 3′-end of a telomere and maintains the telomere during normal cell division (mitosis) (Eisenstein, 2011;Flores et al., 2005). Unusual overexpression of telomerase causes massive extension of telomeric ends and anomalous cell proliferation, which causes cancer. Many studies report that 80-85% of cancerous cells have overexpression of telomerase (Cong, Wright, & Shay, 2002;Lin et al., 2008;Rodriguez-Brenes, & Peskin, 2009). Researchers have hypothesized that the inhibition of telomerase can be an effective approach to return cancer cells to natural apoptosis (Finkel, Serrano, & Blasco, 2007;Sun & Hurley, 2001;Zaug, Podell, & Cech, 2005). The finding that cells switch to alternative mechanisms to maintain telomeric extension is valuable for adopting different strategies for the down regulation of genes (Hu et al., 2012;Ward & Autexier, 2005).
The overhanging G-rich repetitive DNA units at the 3′-end of telomeres can form various tertiary structures, including G-quadruplexes, where the guanine bases stack over each other and are stabilized by the cyclic Hoogsteen type of hydrogen bonding (Burge, Parkinson, Hazel, Todd, & Neidle, 2006;Dai, Carver, Punchihewa, Jones, & Yang, 2007;Luu, Phan, Kuryavyi, Lacroix, & Patel, 2008). These unusual but functionally important DNA structures are also present at other genomic positions, such as promoter oncogenic regions, including c-kit (Fernando et al., 2006), c-MYC (Jain, Grand, Bearss, & Hurley, 2002), bcl-2 (Dexheimer, Sun, & Hurley, 2006), VEGF (Sun, Guo, Rusche & Hurley, 2005), and others. There has been immense research on the biological and physical properties of these noncanonical DNA structures; this research indicates that G-quadruplexes might be a target for designing new anticancer therapeutics (Balasubramanian & Neidle, 2009;Neidle, 2009;Neidle & Parkinson, 2008;. The most promising aspect of the G-quadruplex structure is its topology, which cannot be recognized by the single-stranded RNA component of the telomerase enzyme (Rodriguez et al., 2008). In contrast, the G-quadruplex itself can be a signal of DNA damage and can, therefore, induce apoptosis (Rodriguez et al., 2008). Currently, researchers are more focused on developing agents that can stabilize the secondary G-rich structures of DNA rather than developing agents to inhibit telomerase activity (Dailey et al., 2009;Monchaud & Teulade-Fichou, 2008).
Determining which organic chemical class should be investigated experimentally is a crucial procedure in the initial stage of drug development. The synthesis, derivatization with modification in functional groups, and addition-deletion of excluded volumes for achieving a molecule that is truly specific for the target are indeed tedious, requiring a huge investment of manpower and finances. Therefore, despite much investment in research, no drug has been developed as a potential stabilizing agent of G-quadruplexes. Alternatively, structure-based drug design is an efficient method to screen several potential chemical classes of compounds; after lead optimization, the resulting chemical moiety may prove worthy to the pharmaceutical community (Holt, Buscaglia, Trent, & Chaires 2011;Neidle & Thurston, 2005). Because of the advancement of computing power, in silico drug design has proven to be advantageous in the rational design of novel therapeutic agents; this method is currently in high demand. In the present study, we have adopted the most consistent and efficient methodology of molecular docking to screen a drug library for potential G-quadruplex-stabilizing agents. We have attempted to validate the methodology at regular intervals and to reduce the selection criteria for molecules to produce the most accurate positive result. The docking result was further strengthened theoretically with another promising method, Molecular Dynamics (MD) simulation. In addition, the results of the theoretical studies were validated using low-resolution and high-resolution spectroscopy. Fluorescence and Nuclear Magnetic Resonance (NMR) spectroscopies were employed to show that the hit molecules bind selectively to the G-quadruplex structure. The fundamental nature and various other characteristics of the hypothesized molecules, such as electrostatic properties, solvation free energy, and pharmacokinetics (Absorption, Distribution, Metabolism, Elimination and Toxicity [ADMET]), have been analyzed; this information could be used by various research communities to test and develop these G-quadruplex stabilizing agents.
We search for a promising algorithm that can be adapted successfully to satisfy our goal of finding a G-quadruplex stabilizing agent via in silico docking/virtual screening. All of the available programs have cuttingedge results against protein targets, yet their efficacies against nucleic acids are still debated. Anthony and coworkers have performed docking simulation experiments where the small molecules are targeted to the minor groove of DNA (Anthony et al., 2005). The simulation showed an excellent agreement with experimental results from X-ray crystallography, NMR, and gel-based footprinting. In another type of assessment, it was found that protein-based docking programs such as Glide and GOLD are efficient enough to reproduce the binding modes of the 60 tested RNA complexes (Li et al., 2010). A similar type of validation was also reported recently by a molecular modeling group; they proved that the performance of Glide and GOLD with the simulated structure is very good and produces accurate results in accordance with the structures elucidated by spectroscopic techniques (Srivastava, Chourasia, Kumar, & Sastry, 2011). Moreover, a comparative analysis of different DNA binding drugs for disease targets such as leishmania has also been reported where Glide is used for docking studies with nucleic acid targets (Chauhan, Vidyarthi, & Poddar, 2012). Considering these reports, we have adopted the efficient Glide docking algorithm to screen a drug library for potential G-quadruplex stabilizing agents.

Chemicals and reagents
The DNA duplex and G-quadruplex sequences were purchased from Kric, India. The small molecules, A and B, were purchased from Fisher Scientific Inc (Maybridge, Cornwall, UK). Deuterium oxide was obtained from Cambridge Isotope Laboratories Inc.
G-quadruplex and molecular data-set For this investigation, the NMR structure of human telomeric G-quadruplex TAGGG(TTAGGG) 3 (PDB accession code: 2ld8) (Heddi & Phan, 2011) was chosen. The potassium (K + ) ion was described as an important factor contributing to the structural topology of G-quadruplexes (Lim et al., 2009;Renčiuk et al., 2009). Two K + ions were manually placed between three quartets using the edit option in Maestro. Potential binding sites for G-quadruplex were identified using SiteMap v2.5; the calculations begin with an initial search stage that finds one or more regions called sites, which are suitable for binding to ligand molecules. The molecular data-set containing the small molecules used for virtual screening was obtained from the Maybridge database (HitFinder) of 14,400 molecules. This database was chosen because it has selection of molecules using a clustering algorithm that employs standard Daylight Fingerprints with the Tanimoto (Butina, 1999) similarity index. Clustering with .71 similarities can then be used in NMR for screening against the target (Baurin et al., 2004). Moreover, the entire molecular data-set is in the permissible range of drug-likeness, i.e. the Lipinski rule of five (Lipinski, Lombardo, Dominy, & Feeney, 2001): ClogP ≤ 5, H-bond acceptors ≤ 10, H-bond donors ≤ 5, and mol. wt. ≤ 500.

Docking Studies
All of the molecules were docked using Glide v5.7. The conformational flexibility of the ligand molecule is determined by an extensive conformational search augmented by a heuristic screen that eliminates unsuitable conformations with long-range internal H-bonds. For each core conformation of a ligand, an exhaustive search for possible locations and orientations is performed over the surface of the G-quadruplex. The refined ligand undergoes energy minimization based on precomputed OPLS-AA (Jorgensen, Maxwell, & Tirado-Rives, 1996) van der Waals and electrostatic grids to reduce the large energy and gradient terms resulting from inter-atomic contacts that are too close. A grid with dimensions of 40.10 Å Â 38.63 Å Â 30.56 Å over the G-quadruplex was used to define the binding region and to allow ligands to dock at distances ≤ 20 Å. The defined grid covers the entire topology of the G-quadruplex so that the ligand position and orientation relative to the receptor can be sampled sufficiently; the conformation of the receptor was also fixed during the docking. The standard precision (SP) mode of the software was used for the initial screening of the molecules. Selected top scoring molecules binding to the desired core of the G-quadruplex were selected for the extra precision (XP) mode of docking (Eldridge, Murray, Auton, Paolini, & Mee, 1997). To validate the docking results, we used the docking score of Quarfloxin (CX-3543) (Simonsson et al., 1998), a lead molecule designed and synthesized by Cylene Pharmaceuticals, as a positive control for the docking procedure. Quarfloxin, a pentacyclic system with amino side chains and high selectivity for G-quadruplexes, is already in a phase II clinical trial. Other classes of compounds with reported interaction with G-quadruplexes were also used to perform the control test of docking studies (Zhen, Xi, Chattopadhyaya, & Hecht, 2011). The Glide score is based on the ChemScore scoring function coupled with the ligand-receptor molecular mechanics interaction energy and ligand strain energy; thus, provides a docked pose. Based on the binding site of these drugs and their docking score, we then validated our docking study of small molecules from the Maybridge database. As part of the docking result analysis, we also used a 1000 ligand-like decoy set from Schrödinger (www.schrodinger.com/productpage/14/5/74/) to evaluate the results using a receiver operating characteristic curve (ROC) curve approach (Triballeau, Acher, Brabet, Pin, & Bertrand, 2005). The characteristic area under the curve obtained from ROC can be helpful in predicting whether the obtained hits are true positives.

Molecular dynamics (MD) simulation
Dynamic simulation of the drug-bound and drug-free states of the G-quadruplex was performed using Desmond Molecular Dynamics System v2.2 with OPLS-AA force fields (Jorgensen et al., 1996). The OPLS force field generates accurate charges and force field parameters for ligands. This force field was adopted with an objective to keep the same set of protonation states and charge parameters for ligand pose refinement and conformational sampling during docking; thus, the same parameter set will also be used in the simulation studies. The solute system containing the drug-bound and drug-free states of the G-quadruplex was solvated using the TIP3P water model (Jorgensen, Chandrasekhar, Madura, Impey, & Klein, 1983), with orthorhombic (α=β=γ=90°) boundary conditions extending a distance of 10 Å from any solute atom in all directions. The charge of the total system was neutralized by adding the appropriate number of Na + counter ions. The number of explicit water molecules ranged from 3891 to 4008. After the set up of the system, an initial minimization was performed with a convergence threshold of 1.0 kcal.mol À1 Å À1 to allow all of the atoms to adjust to the system environment. An MD simulation was performed for each system using an isothermal-isobaric ensemble (NPT) with a relaxation process. The relaxation process consisted of minimization with and without solute restraints and simulation in an canonical ensemble using a Berendsen thermostat with a temperature of 10 K and restraints on solute heavy atoms for 12 ps. In contrast, simulation in the NPT ensemble was performed using a Berendsen thermostat and a Berendsen barostat with a temperature of 10 K and a pressure of 1 atm for 12 ps. In the next stage of relaxation, the system was simulated for 24 ps with Berendsen NPT with a temperature of 300 K and a pressure of 1 atm with and without restraints on solute heavy atoms. The particle mesh Ewald (Darden, York, & Pedersen, 1993) method was used with a tolerance of 1e À09 as a part of the Coulomb method. A cutoff of 14 Å was used for calculating the nonbonded solvent-solvent and solute-solvent interactions. The pressure and temperature of the system was maintained at 1.013 atm (isotropic) and 300 K, respectively. The SHAKE algorithm (Krutler, Gunsteren, & Hnenberger, 2001) was employed with an integration time step of 2 fs. A recording interval of 1.2 and 4.8 ps for the energy and trajectory frames, respectively, was set. The MD simulation for the entire system was continued for up to 5 ns, and the trajectory analysis was performed using the Simulation Event Analysis module of Desmond (Scheme 2).

Adaptive Poisson-Boltzmann Solver (APBS) iso-surface calculations
The Adaptive Poisson-Boltzmann Solver (APBS) was used to analyze the electrostatic characteristics of the macromolecular complex. The standard all-atom charges generated by the OPLS force field were used for this purpose. The Poisson-Boltzmann equation (PBE), the most widespread model for evaluating electrostatic properties, was adopted from the reference of Baker and coworkers (Baker, Sept, Joseph, Holst, & Cammon, 2001). The essence of the model is that it relates the electrostatic potential to the dielectric properties and the distribution of the atomic partial charges of the biomolecular system. The various computations required to evaluate the electrostatic properties were performed using a plugin option of PyMOLv1.2r3pre, and the resulting electrostatic contour visualizations were collected at the appropriate positive and negative iso-surface values, optimized for visualization.
Free energy perturbation (FEP) and concept of free energy The main objective of performing the MD simulations for an ensemble system is to explore the macroscopic properties of that system based on its microscopic properties. The macroscopic properties include various thermodynamic terms, such as temperature, pressure, number of particles in the system, and the free energy of binding of a particular compound; in contrast, the atomic positions (r) and the momenta (p) are the microscopic properties of a system, which are defined as the coordinates of the multi-dimensional space called the phase space. The concept of free energy calculations was developed in the early nineties (Kirkwood, 1935;Zwanzig, 1954) but was not popular because of the complexity. The first successful treatment of free energy calculations was for the proton transfer in the lysozyme (Warshel & Levitt, 1976). These calculations gained popularity with significant works by several scientists (Kollman, 1993;McCarrick & Kollman, 1999) who integrated thermodynamic parameters with the free energy calculations. This free energy calculation is mainly based on Zwanzig's equation, which gives the free energy difference between the two states A and B.
where, β = 1/kT (k is Boltzmann constant and T is temperature), […] A denotes an MD or Monte Carlo generated ensemble average of state A, and ΔV = V A ÀV B . V A and V B are the potential energies of A and B states. For absolute solvation free energy expressions, A can be defined as the solute molecule in the gas phase and B as the solute molecule in the solution phase. The free energy calculations were performed in Desmond. From Figure 1, the solvation free energy for Figure 1. Schematic representation of the free energy change in solute-solvent interaction in different phases. Reactant is referred to ligand molecule and since there is no differentiation in the formed product except the molecular conformation, nothing is referred to the products here. ligand molecules can be calculated by annihilating the reactant (ligand) in solvent and in vacuum.
The main advantage of Desmond is that it requires only one step to calculate the solvation free energy of the solute; the reactant in vacuum is recovered from the solvent phase by rejecting the interactions of reactant with any other atoms in system (Shivakumar et al., 2010).
The solvation free energy calculation plays a major role in predicting the binding affinity of the compound and is a powerful tool in the lead optimization step for drug discovery. Implementing a FEP is a multi-step method, which involves the use of a group of intermediate potential energy functions defined between the paths from state A to state B (Brandsdal et al., 2003). In these functions, λ is the coupling factor whose value ranges from 0 to 1, and m = (1 to n) is the number of points; each point possesses an individual potential energy function corresponding to the particular λ value.
A new process of multiple step free energy calculations was proposed, which is called "λ-dynamics." In this process, the λ value is treated as a dynamic variable whose value ranges from 1 to n (Kong & Brooks, 1996).
The main concept involved in this method is that different ligands compete for the same receptor in the body; this competition can be modeled by the simultaneous searching of the various ligands with a varying value of λ (Elber & Karplus, 1990). The use of various intermediate states increases the accuracy of the free energy calculations. A more accurate method for the calculation of free energies was proposed by Bennett (1976) and referred to as the Bennett Acceptance Ratio. The Free energy perturbation (FEP) calculations were performed in Desmond using this method in which a number of simulations were run in parallel for a single FEP calculation.
For the solvation free energy calculation in Desmond, the total free energy by the FEP workflow was employed: the small molecules were annihilated separately in the solvent system with the TIP3P water models, the simulation was run for up to 2 ns, and λ was given a value of 11. A correction value for the energy, which represents the misplaced long-range dispersion energy because of the cutoff of the van der Waal potential energy, was also obtained using FEP calculations.

Fluorescence studies
Fluorescence emission spectroscopy was used to analyze the binding behavior of both molecules (A and B) with a G-quadruplex, 5′-TAGGG(TTAGGG) 3 -3′. A G-rich DNA duplex was used as a control in this study with the following sequence: 5 0 À GCGCATGCTACGCG À 3 0 3 0 À CGCGTACGATGCGC À 5 0 . Emission spectra were recorded at room temperature using a Hitachi spectrometer (F-700 FL spectrophotometer) with 1 cm path-length quartz cuvette. Each emission spectrum was recorded from 300 to 500 nm with an excitation wavelength of 260 nm. The excitation and emission slits were set to 5 nm each. The emission spectra were recorded by titrating the small molecules with concentrations of up to 50 μM against 10 μM of the G-quadruplex and the duplex DNA. All fluorescence experiments were performed in 10 mM phosphate buffer, 10 mM EDTA, and 100 mM KCl (pH 7.0). Because of solubility problems with molecule A, a small amount of NaOH was added to dissolve it in at the same buffer.
WaterLOGSY NMR experiments All spectra were recorded at 298 K with a Bruker AVANCE III 500 MHz NMR spectrometer equipped with a 5 mm smart probe. For each sample, a reference spectrum and a 1D WaterLOGSY spectrum were recorded in 10 mM phosphate buffer, 10 mM EDTA, and 100 mM KCl (pH 7.0) containing 10% D 2 O. The first water-selective 180°pulse was 25 ms long. A weak rectangular PFG was applied for the entire mixing time (1.5 s). A short gradient recovery time of 1 ms was introduced at the end of the mixing time before the detection pulse. The two water-selective 180°square pulses of the double spinecho sequence were 2.6 ms long. The gradient recovery time was 0.2 s. The data were collected with a sweep width of 10 ppm, an acquisition time of 1.5 s and a relaxation delay of 2.0 s. Prior to Fourier transformation, the data were multiplied with an exponential function with a line broadening of 2.0 Hz. All chemical shifts were referenced using (2,2-Dimethyl-2-silapentane-5-sulfonate sodium salt) as an internal standard.

Results and discussion Docking Studies
Two types of binding sites were found for the G-quadruplex from SiteMap analysis (Figure 2). By structure, these two binding sites can be classified as: (a) a loop region, which formed small cavities and (b) a co-facial/ end stacking region at the top of the macromolecule (Scheme 1). The loop regions are highly rich in electron density and are favorable for binding with small molecular fragments. These loop regions provide good binding affinity for hydrogen bond acceptors. The other binding site, i.e. the co-facial end stacking region containing the stacked guanine residues (G-face) (forming the first quartet), is enriched with π-electron clouds from all four directions (Scheme 1). This site is important for the design of novel G-quadruplex stabilizing agents because it is the most suitable region for binding of small molecules, peptides or DNA-RNA aptamers (Murat, Singh, & Defrancq, 2011), which can add favorable entropy to the system, resulting in its stabilization. The region was explored as the most suitable region when we docked the previously reported stabilizing agents with the G-quadruplex ( Figure 3); hence, we excluded the investigation with the loop binders. This initial docking study not only enables us to explore the favorable binding site for interaction but also provides a means to validate our docking studies (Scheme 1). We screened a set of 14,400 molecules that were available in the Maybridge database using the SP mode of Glide.
The docked poses obtained for the molecules from this study were superimposed with the docked poses obtained for the test set of molecules. The selection criteria for the molecules to be processed in the XP mode of docking were that the Glide SP score was comparable to that of the test set with a rmsd fit of up to 2 Å for their binding pose. A total of 42 molecules of comparable Glide score were selected; these molecules had a docked pose in the co-facial region only and did not have any unusual binding poses in the loop region. Although the basic docking algorithm utilized by Glide is the same for both SP and XP modes, the alternative sampling of ligand conformation with ligand-receptor complementarity and scoring function is more extensive in the XP mode. Hence, with the XP mode of docking, the false positives can be excluded. The top Glide score obtained from the XP mode of docking was À5.65 for Quarfloxin (the docked pose is represented in Figure 3). Two molecules, A and B, showed satisfactory results in agreement with the score and conformation of Quarfloxin. The two molecules, N′1-imino(2-pyridyl)methyl-3,4,5-trimethoxybenzene-1-carbohydrazide (molecule A) and ( Figure 4. The competence of the docking result was also assessed with a set of ligand decoys containing 1000 molecules. These decoys are nonspecific binders and can bind to any part of the receptor macromolecule. This data-set is known to consist of inactive molecules; they will be treated as false positive even if they bind to the receptor site. Thus, the docking efficiency and accuracy of the 42 selected molecules can be compared with respect to their binding pose using a receiver operating characteristic curve (ROC) ( Figure 5). The ROC value for the plot, "sensitivity" vs. "1-specificity", was .865; interestingly, the area under the fitted curve was found to be .850. This result dem- onstrates that the docking was accurate. The high ROC value (close to 1) indicates that the 42 selected molecules in our study selectively bind to the G-core of the quadruplexes (Table S1).

Structural stability of drug-bound G-quadruplexes
The high docking score for two molecules, A and B, suggested that they may have a high affinity for G-quadruplexes. To generate the additional theoretical information required to design a new therapeutic agent, a 5 ns MD simulation was performed with both ligand molecules bound to the G-quadruplex. The simulation fetch average conformations of the drug-like molecule bound to the G-quadruplex with very few fluctuations in rmsd    (Table S2). A free-state G-quadruplex was also simulated with the same time scale and used as a control and reference for the simulation results of the drug-bound state. Figure 6 (A-C) shows the rmsd plot for all of the models compared with the initial structure and considering various atoms under investigation. In the ligand-bound state, the fluctuations, based on the rmsd, of the backbone are higher than those for the nucleotide bases. The greater variation in the rmsd of the free-state G-quadruplex than that of the ligand-bound states indicates the absence of steric strain by ligands. Ensemble structures of various MD frames are aligned and shown in Figure 7. The rmsd deviation of molecules A and B bound to G-quadruplexes appears to be almost comparable for the bases and for the backbone; this finding demonstrates that both ligands behave as stabilizing agents with prominent contacts and are potent in controlling the fluctuations of G-quadruplex structures (Table S2). An initial jump in the rmsd scale during first few ps of simulation compared with that of the starting frame corresponds to the relaxation of the model. The stability of the trajectories of the drug-bound states can be observed from 1 to 5 ns, with minimal fluctuations ranging from 0.5 to 1.5 Å. The rmsd of both drug-bound states are limited to this deviation and indicate that the average structures of the bound states are approximately convergent; thus, the MD simulation was terminated after 5 ns. In contrast, the rmsd plot for the free-state G-quadruplex does not appear to have stabilized in 5 ns ( Figure 6). Although achieving stability in the free-state model was not the subject of investigation, it was used as a control for the study of A and B in their bound states. A plot of the allatom rmsd vs. time also explains the similar behaviors ( Figure 6(C)). The standard deviations for the all atoms rmsd of A and B were $0.487 and $0.427 Å, respectively, whereas that of the free-state G-quadruplex was $0.697 Å .

Coulombic and van der Waal interaction energies
The van der Waal and Coulomb energies are two important contributors to the stabilization of any biomolecule. Table 1 describes the van der Waal and Coulomb interaction energies (kcal.mol À1 ) of G-quadruplex with selfatoms, solvent atoms (water models), and solute atoms of ligands and ions. The Coulomb energy is the major contributor for the stabilization of the drug-bound state of G-quadruplex (Table 1). The Coulomb energy for the model favors binding to the ligand compared with the free-state of the model. The strain of ligand molecules makes the aromatic nuclei of guanine bases stack over each other, which can be judged by the average torsion angle deviations (Table S2). Furthermore, the protonation by the ligand compensates for the electronegativity of O6 by restricting its fluctuations. When shared, the electron orbitals of O6 might become limited to a certain number of degrees of freedom, which ultimately limits the deviation in geometry of the guanine core. This type of effect was shown by comparing the rmsd of O6 specifically over time in all three models, as shown in Figure 8. The central N atom of molecule A is able to   Figure S1 in supporting information). The same difference can be observed in the corresponding energy contributions shown in Table 1. In molecule A, there are three methoxy groups and one of these methyl groups interacts with the π cloud of the guanine ring. This type of σ-π contact is very favorable and is another reason for stable binding of molecule A-G-quadruplex. A diagrammatic explanation is given in the supplementary data ( Figure S2 in the supporting information). The Coulomb and van der Waal contributions to the stacking energy are almost comparable for molecules A and B. As in the free state, there is no ligand interaction with the G-quadruplex; the contribution from Coulomb interaction is almost zero. Therefore, the total energy comes from the van der Waal interactions between two K + ions and the G-core. The majority of the effects of the continuum system of the TIP3P water model are from the Coulomb interactions rather than the van der Waal interactions. This interaction is more pronounced in the absence of a ligand. A similar type of explanation is given using a Poisson equation describing the multi-pole representations of the electronic distribution in Figure 9.

Hydrogen bonding analysis
Hoogsteen type H-bonding is the key factor behind the topology of the G-quadruplex. The average inter-atomic distances between N1H-O6 and N2H-N7 are shown in Table 2. The hydrophobic core of the first quartet is the most suitable region for the binding of stabilizing agents. Searching for structural deviations in the inter-atomic distances of polar atoms can provide valuable information regarding the interaction of DNA with small molecules. One terminal of molecule B was buried sideways between the first and second quartet, binding it to the model with a strong H-bond, and the other terminal was above the structure ( Figure S3 in the supporting information). The Hoogsteen type H-bonding plays a vital role in the stabilization of guanine structures. The planarity and uniformity of the H-bonds involved in Hoogsteen G-pairing is well controlled by the position of the O6 atom. A plot of rmsd (O6) vs. time for the drug-bound states of the G-quadruplex showed that the displacement of the O6 atom from the plane of the Hoogsteen base pairing is negligible during the 5 ns of simulation ( Figure 8). For the ligand-receptor complex, we found that the Hoogsteen H-bonding in the G-pairing and H-bonding between the ligand and receptor play a key role in stabilizing the complex; πÀπ stacking interactions also enable the correct positioning of the ligands over the G-quadruplex for the entire time scale. As indicated in Figure 4, the atoms of ligands that are emphasized with the arrow H-bond to the macromolecule. The two hydrogen atoms that have a central protonated N atom (H-bond donor) in molecule A make hydrogen bonds with the O6 of guanines (G5 and G23) ( Figure S1 in the supporting information). In contrast, one terminal N atom, a H-bond acceptor, from a cyanide (-CN) group of molecule B makes a H-bond with the amine group (-N(2)H) of a guanine residue (G4) in the second quartet ( Figure S1 in the supporting information) (below the facial quartet).
The average H-bond distances computed for molecules A and B were $2.54 and $2.45 Å, respectively, which emphasizes that H-bonds are the main contributor to the binding of both ligands to the G-quadruplex. The strength of the H-bonding for the H-bond acceptor in molecule B is much greater than that of the H-bond donor of molecule A. This finding can be seen in a histogram plot (Figure 10), which counts the frequency of H-bond distances for both molecules. The Hbonding pattern provides a steric restriction on the rotational freedom of the C-C bond in the ligands. This restriction in the rotational degrees of freedom reinforces the stacking of the aromatic moiety of the ligand over the G-core. The H-bonding arrangement of molecule B is biased toward the side of the molecule in which the geometry is most affected by the conformational constraints of the ligand-DNA complex. This bias comprises the flexibility to its other half; the dynamicity is shown in a plot of distance vs. time (Supplementary data Figure S3). Examining the overall structure indicates that a small molecule should be rationally designed with a N atom that could be utilized as a donor (molecule A) or acceptor (molecule B) to form a contact with the G-core; doing so can strengthen strong binding in the system.

Electrostatic Iso-surface by APBS
A standard practice in molecular biophysics is the evaluation of the electrostatic properties of a ligand-macromolecule complex. The APBS, developed by the group of Nathan Baker, is a routine using the PBE to illustrate the electrostatic attributes of a biomolecular complex; thus, the molecular electrostatic surface can be understood and used to develop a good stabilizing agent for G-quadruplexes. Figure 9(A) shows the electrostatic isosurface contours of the G-quadruplex model with positive and negative iso-surfaces at ± 4 kT/e. If we carefully analyze the G-quadruplex model, the macromolecule can be assumed to be a system with highly anionic characteristics. For the ease of comparison, a top view of the G-quadruplex model is shown with the phosphorous atoms in sphere representation (Figure 9(B)), indicating the relative position of the phosphate groups, sugars, and bases in the quadruplex topology. In the G-quadruplex structure, the electronegative charge is diffused in two ways: (1) the negative charges of the backbone phosphates surround the nucleobases and sugar moieties and (2) the lone pair electrons of guanine O6 atoms are directed toward the center, increasing the electron density in the G-core. Because of the dense cloud caused by negative charge diffusion, the iso-surface contours of the G-quadruplex form a large dome shape. The π-electrons of the purine rings also add to the overall electronegativity of the system. When we compared the electrostatic properties of the G-quadruplex system bound to the drug, as shown in Figure 9(C) and (D), we found a reduced overall charge in the system in the hydrophobic region. Because of this reduction, a G-quadruplex can only bind to (and thus recognize) a molecule if it has sufficient positive electrostatic fields. In molecule A, the methyl protons also interact with the aromatic system of the preceding guanine base, which tightens the binding of molecule A with its host. The electrostatic iso-surfaces of the G-quadruplex bound to both molecules A and B were obtained at ± 4 kT/e. The electrostatic iso-surfaces of both ligand systems are shown in Figure 11, where the iso-surfaces were achieved at ± 2 kT/e, to illustrate the broad range of electropositive values. The partial atomic charges of individual molecules were assigned according to the OPLS force-field and are shown in Figure 4.

Solvation free energy
The solvation energies of molecules A and B were computed in pure solvent. The annihilation of the compounds resulted in a total solvent energy of À48.7 ± 0.3 kcal. mol À1 and À74.5 ± 0.6 kcal.mol À1 for molecules A and B, respectively. These values are the lowest terms of the 11 λ windows that were separately annihilated to find the various possible low energy conformations in which the ligands were expected to bind G-quadruplexes. Because the ligands contained no intra-perturbed groups for the solvation free energy calculation, the interaction of both of the small molecules with the solvent must be favorable. In molecule A, the central charged amine group is bound to an aromatic system (pyridyl group) from one side, which hinders the solvation shell of the amino group. Similarly for molecule B, the two charged central amino groups are affected by a benzylidene group that utilizes the lone pair availability of the N atom and thus increases the solvent accessibility. This effect can be understood in terms of the solvation free energy difference between molecules A and B. Graphical plots  To the best of our knowledge, there is no report of a study emphasizing the calculation of the standard molecular solvation free energy of Gquadruplex stabilizing agents. Without this experimental value, it is very difficult to validate our theoretical data. When analyzing the plots for both ligands for the solvation free energy with the same time scale, molecule B shows a relatively higher free energy. This finding Figure 11. Iso-surface contours for molecule A and B calculated at ± 2kT/e. Positive and negative isosurfaces were represented by blue and red colors, respectively.
indicates that molecule B is more solvated than molecule A. These data are based on the larger solvent accessible surface area (SASA) of B and attributed to the higher degree of freedom in the state bound to the G-quadruplex. The solvation free energy calculation is also accompanied by a correction factor, which is typically known as a long-range dispersion energy correction. The correction values were À3.0 and À4.2 for molecules A and B, respectively. The details of the solvation free energy calculations are given in Table 3. Few of the other relevant properties computed from QikProp can be correlated with the discussion of solvation free energy in       Table 5. Both compounds are Central Nervous System (CNS) -inactive; therefore, it can be hypothesized that the molecules have no central nervous system activity. All of the other predicted properties described in Table 4 appear to be within the permissible range.
Binding of the G-quadruplex with small molecules measured using fluorescence spectroscopy To validate the theoretical study, we used the intrinsic fluorescence of the G-quadruplex to monitor its binding with molecules A and B. A significant hypsochromic shift or blue shift in the fluorescence emission provides strong evidence for whether a drug/ligand molecule binds to the G-quadruplex. The fluorescence experiment was performed by titrating molecules A and B (0-50 μM) with the G-quadruplex, 5′-TAGGG (TTAGGG) 3 (10 μM). The characteristic emission of the spectrum is shown in Figure 13. A total of 10-18 nm hypsochromic shifts in the G-quadruplex emission spectra were observed when either molecule A or B was bound to the G-quadruplex (Figure 13(A) and (C)). The changes in the fluorescence emission maxima of the Gquadruplex as a function of the concentration of the ligands, molecules A and B, yields equilibrium dissociation constants (K d ) of 137 ± 1.1 and 31 ± 1.2 μM, respectively ( Figure 13(B) and (D)). Molecule B is a better binder to the G-quadruplex than molecule A. A control experiment was performed for a GC-rich DNA duplex, 5 0 À GCGCATGCTACGCG À 3 0 3 0 À CGCGTACGATGCGC À 5 0 in the presence of molecules A and B. Interestingly, no blue shifts of the emission spectra of the DNA duplex occurred with the addition of the ligands (Figures 13(E) and (F)), but a significant increase in the intensity in the fluorescence spectra upon the addition of the ligands was observed. To determine the cause of the unexpected increase in fluorescence intensity, we measured the fluorescence spectra of only the ligands, A and B, without any interference of the host, the G-quadruplex ( Figure S4 in the supporting information). We observed an increase in fluorescence intensity in the emission spectra at $340 nm for both molecules when they were excited at 260 nm ( Figure S4 in the supporting information). This finding explains why the G-quadruplex and DNA duplex showed increases in fluorescence intensity upon the addition of ligands A and B. However, experimental evidence such as blue shift of G-quadruplex upon addition to the molecules A and B clearly showed that molecules A and B bind much more selectively to the G-quadruplex structures than to the DNA duplex.
Water ligand observed via gradient spectroscopy (WaterLOGSY) NMR experiments Water ligand observed via gradient spectroscopy (Water-LOGSY) NMR experiments allow the detection of binding of small molecules to protein targets with dissociation constants in the low micromolar to millimolar range (Dalvit et al., 2000;Dalvit, Foqliatto, Stewart, Veronesi, & Stockman, 2001). Compared with a saturation transfer difference NMR experiment (Bhunia, Saravanan, Mohanram, Mangoni, & Bhattacharjya, 2011;Mayer & Meyer, 1999), which is also an important technique for measuring the binding of small molecules to receptor proteins, WaterLOGSY is more effective for studying the interactions of nucleic acids and small molecules (Bhunia, Bhattacharjya, & Chatterjee, 2012;Lepre, Moore, & Peng, 2004). The WaterLOGSY signal is produced by the selective transfer of bulk water magnetization to the ligand in solution via the macromolecule-ligand complex. In the resulting spectra, the ligands that bind generally have an opposite sign to those of compounds that do not bind. Figure14 shows the one-dimensional WaterLOGSY NMR spectra of the small molecules, A and B, in the presence of the G-quadruplex. Interestingly, both molecules bind to the G-quadruplex, as indicated by the positive signals. The aromatic ring protons and the methyl protons of molecule A show positive signals in the presence of the G-quadruplex. The simulation data suggests that the three methoxy groups attached to the benzene moiety are in close proximity to the G5 and G23 groups of the G-quadruplex. The aromatic ring protons of benzene and the pyridyl groups of molecule A severely overlap at 7.35 and 7.52 ppm (Figure14B). However, the MD simulation data indicate that the aromatic rings, benzene and pyridyl, are stacked to G11, G5, and G17, respectively. The putative binding site and the MD simu-lation data correlate well with the WaterLOGSY data ( Figure 14(C)). In contrast, the intensity of the peaks of molecule B is significantly stronger than that of molecule A, when binding to the G-quadruplex. Interestingly, all of the protons, except for those in two ethylene groups of molecule B, show positive WaterLOGSY signals. This result demonstrates that the ethylene groups (3.12 and 3.82 ppm) do not interact with the G-quadruplex as indicated by the negative signals (Figure 14(E)). This finding is not surprising because the MD simulation data clearly show that one terminal of molecule B is quite flexible and that the ethylene groups point toward the solution (Figure 14(F)). Both N-Me groups (3.67 and 3.38 ppm) and both aromatic ring protons (6.90 and 7.83 ppm) and the flexible linkers in the center show strong Water-LOGSY signals. In the simulated model, G17 and G23 are in close proximity to one terminal of molecule B, whereas G5, G10, and G23 are very close to the other terminal of the G-quadruplex (Figure 14(F)). The flexible linker, attached to both aromatic ring protons of molecule B, stacked to the G-faces, such as G11, G17, and G5 of G-quadruplex. As a control experiment, we measured the WaterLOGSY signals of molecules A and B in the absence and presence of GC-rich duplex DNA (Figure S5 in the supporting information). None of the molecules bound to duplex DNA structures. The WaterLOGSY NMR experiments in conjunction with the MD simulation characterize the binding epitope of both molecules, which bind to the G-quadruplex at an atomic resolution.
A final hypothesis: requirements for the rational design of G-quadruplex stabilizing agents: From previous literature and this study, we elaborate some of the key features, which can be implemented for the future rational design of a G-quadruplex stabilizing agent. In designing a suitable stabilizing agent for guanine-rich secondary structures, the primary concern is to understand the structural criteria needed to complement the G-quadruplex. Some of the criteria for a rational design of novel G-quadruplex stabilizing agents are as follows: (1) A molecule with enhanced aromaticity is essential to effectively interact with G-quadruplexes. π-π stacking interactions are key requirements for any molecule to be attached to the facial guanine quartets.
(2) H-bond acceptability is the main parameter for a ligand to bridge with the loop surrounding the Gcore. In molecule B, nitrogen can act as a Lewis base and donate its lone pair of electrons to form a hydrogen bond, which helps place molecule B above the quadruplex structure. This property is important because it establishes a criterion to develop stabilizing agents irrespective of the facial topology of G-quadruplex. (3) In contrast, the donation of a H-bond to the quartets by a ligand molecule is a prerequisite for the facial G-quadruplex stabilization. Ligand molecules with this characteristic will enable the utilization of the lone pairs available from the O6 of guanines. This effect is more pronounced when the H-bond donor moiety is found in the central region of the ligand molecule, as in molecule A. Apart from this, the protonation to this moiety provides a positive potential directed toward the central hydrophobic region. (4) The methyl-π type interaction also helps to stabilize the ligand and G-quartet assembly. Ligands with a methyl moiety localized on the top of the central part of the G-quadruplex try to stabilize the host-guest assembly with methyl-π type dipole-induced dipole interactions, as is the case for molecule A. This type of nonpolar interaction from the mixing of orbitals of the methyl-proton with the π electron in the purine ring of guanine is of higher order. This interaction can provide a very stable and strong binding affinity with a minimal margin of excluded volumes to any organic core molecule. (5) The radius of gyration of the molecule is another important criterion. The average distance between the guanine residues of a quartet from all types of G-quadruplex structures is within a range of 15-17 Å. If the radius of gyration of the molecule falls in this range, it can successfully find a conformational orientation over the hydrophobic core of the G-quadruplex. Additionally, the middle portion of the molecule should have rotational degrees of rotation with respect to the C-C bond to be able to adjust to the planarity of the guanine face. The radius of gyration of a molecule is related to the dimension of the central core and not of the functional groups attached to the core. (6) Electrostatic interaction has a more significant role than any other physical descriptor. In our study, we found the G-quadruplex structure to be highly anionic with abundant electronic density surrounding the structure from all directions. The incoming ligand must have sufficient cationic characteristics. In other words, the incoming ligand should have adequate electron-withdrawing effects to stay and bind over the planarity of the quadruplex.

Conclusion
In this study, we revealed that the two molecules, A and B, are pharmacologically potent in arresting G-quadruplex structures in DNA. The data from MD simulations indicate that these screened analogs are capable of binding the quadruplex within a reasonable time. Fluorescence experiments were performed by titrating molecules A and B (0-50 μM) with the quartet DNA of sequence, 5′-TAGGG (TTAGGG) 3 (10 μM). A large blue shift was found in the emission spectra of the G-quadruplex (the fluorescence emission of the G-quartet is at 337 nm) with stepwise increase in the ligand molecule concentration over the concentration of the G-quadruplex. Molecules A and B give emission blue shifts of 10-18 nm, respectively, demonstrating that both molecules bind to the Gface of the quadruplex. A one-dimensional WaterLOGSY experiment further proved at atomic resolution that both molecules, A and B, bind to the G-quadruplex structure; the structural information of the bound complex is well correlated with the MDs simulation. Our experimental investigation is underway with these molecules. The framework of these two molecules can be the platform to design novel anticancerous therapeutics in the near future.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.10.1080/07391102.2012.742246.