Skeleton selectivity in complexation of chelerythrine and chelerythrine-like natural plant alkaloids with the G-quadruplex formed at the promoter of c-MYC oncogene: in silico exploration

Abstract

‘Pu27’ is a G-rich 27 nucleotide sequence in the non-coding strand of the ‘nuclear hypersensitive element (NHEIII₁)’, present in the promoter region of the c-MYC oncogene. Pu27 attains a stable secondary structure (G-quadruplex) and regulates the expression of c-MYC, thus stabilization of Pu27 using small molecules has become a novel strategy for developing cancer therapeutics. Recently the plant alkaloid chelerythrine has been explored as a probable anti-cancer agent with its novel mechanism of action i.e. stabilization of the G-quadruplex formed by Pu27. In the present study we have determined the binding characteristics of chelerythrine and three other chelerythrine-like plant alkaloids with the G-quadruplex formed by Pu27. Binding conformations of each ligand are obtained by molecular docking. Furthermore, the structural integrity of the built model and the conformational behavior of Pu27 upon binding of the ligands are established through all atom molecular dynamics simulations in an explicit solvent for 50 ns. Binding free energies of the ligands are estimated over the last 2 ns simulation run. We cultivated the gravity and influence of the position of methoxy group on the chelerythrine skeleton, tuning the binding pattern and energetics over its host, Pu27. Based on the results we categorized the ligands according to their better binding efficacy. Chelerythrine and the 12 methoxy analog of it showed the highest affinity towards complex formation. The in depth binding behaviors of these molecules were extracted from the conformational hyperspace of complexes reflected in their principal component analysis. We report the atomistic details of the binding patterns of chelerythrine and chelerythrine like molecules towards the G-quadruplex formed in Pu27. Chelerythrine and 12-methoxy chelerythrine show similar binding patterns, stacking over the 5′ end of the sequence, methoxy substitution over the 12^th position aids in the binding activity by enhancing the van der Waals energy contribution. A change in the position of methoxy groups over ring-A of chelerythrine alters the ligand binding.

---

1 Introduction

Expression of the c-MYC oncogene has a pivotal role in regulation of cell growth. Its over expression with various mutations is observed in numerous human cancers like breast cancer, colon cancer and many more. c-MYC over expression is closely associated with cell proliferation, it inhibits the differentiation of cells thus causing malignancy.^1–5 It has been studied that the inhibition of expression of c-MYC in cancer cells reinstates the differentiation of cells. Thus transcriptional regulation of c-MYC oncogene is an anticipated mode for anti cancer therapy. ‘Nuclease hypersensitivity element (NHE) III₁’ majorly regulates the transcription of c-MYC through structural transitions of ‘Pu27’;^6–11 Pu27 is a 27-bp long purine rich non coding sequence of NHEIII₁ which is located at the upstream site of P1 promoter of c-MYC. Siddiqui-Jain et al. have analyzed the structural characteristics of Pu27, they have demonstrated the presence of two intra-molecular G-quadruplex structures out of which only one is biological relevant.¹² Also, they have found that the transcription of c-MYC can be suppressed using G-quadruplex stabilizing agents like TMPyP4. Thus ligand mediated stabilization of the G-quadruplex structure of Pu27 is under extensive exploration by means of chemical, biophysical, biological and in silico approaches.^13–18 As such, stabilization of G-quadruplex structures found in the promoter regions of various oncogenes and telomeric region has been considered as an effective measure in anti cancer drug development,¹⁹ our research group has been delving into the structural countenances of G-quadruplex–ligand interactions with the same purpose.^20,21 Also, we have offered the building mechanism behind the formation of G-quadruplex structure.²² Various naturally occurring plant elements and their derivatives are found to be efficient in stabilizing G-quadruplex structures.^{16,18,23–31} Chelerythrine (CHE) is a plant alkaloid with an anti cancer potential; different mechanism of actions are distinguished for its anti cancer activity such as inhibition of protein kinase C,³² Ghosh et al. have established a novel mechanism where CHE induces the secondary structure of G-quadruplex in the telomeric region of human genome and further causes inhibition of telomerase.³³ Moreover they have verified another mechanism of action where CHE binds to G-quadruplex structure of Pu27 and enhances its stability thus down-regulates the expression of c-MYC oncogene.¹⁵ EtBr displacement assay indicated that CHE may interact with Pu27 by π–π stacking interactions over the end regions and further increases the stability of overall secondary structure of Pu27. Till date, the atomistic details of CHE binding to Pu27 is not been explored. In this regard, to unravel the importance of structural architecture of the receptor, Pu27 as a novel target for cancer therapeutics and CHE as its arresting and stabilizing agent, we explored the interaction behavior of CHE with the G-quadruplex structure of Pu27 at the atomic level by means of molecular dynamics simulations in explicit solvent. In addition, we mined out three other plant alkaloids (as shown in Fig. 1) possessing structural similarity with CHE, may have potential to bind Pu27 selectively. These molecules are examined side by side with CHE to decipher the structural as well as functional roles of the methoxy groups at different sites of the ligand. Our key observations could define the pharmacophoric features necessary in stabilizing Pu27. This may assist in the design and optimization of lead compounds targeting Pu27 structure which confers a novel approach towards anticancer therapy.

Fig. 1 2D structures of ligands selected for simulation analysis: (A) chelerythrine (CHE), and other chelerythrine like plant alkaloids named as (B) ONE (12-methoxy CHE), (C) TWO, (D) THREE (nitidine). Respective database IDs for each chelerythrine like compound are also provided. (E) Cartoon representation of model built for Pu27 native sequence.

2 Materials and methods

2.1. Building an initial structure of G-quadruplex for native sequence of Pu27

Previous NMR studies have determined the structural pattern of G-quadruplex formed in the promoter region of c-MYC oncogene. It forms the stable parallel-stranded unimolecular quadruplex structure, all guanines involved in G-tetrad formation are in anti conformation and tetrads are interconnected by propeller loops.^34–36 Pu27 has five guanine tracts in its sequence out of which first guanine tract does not participate in the formation of G-tetrad of respective quadruplex structure,¹² however, it takes multiple secondary quadruplex structures thus gives indistinguishable broadened imino peaks in its NMR spectrum¹⁴ (PDB ID 2A5P). The NMR solution structure of truncated form of Pu27 i.e. Pu24I has four guanine tracts eliminating the first guanine tract of Pu27, these four guanine tracts form three stacked guanine tetrads with the parallel arrangement and all stacked guanines attain anti conformation. Though the first guanine tract is not necessary for structural features of respective G-quadruplex but it may influence the binding of stabilizing agents thus it is essential to consider the entire sequence for ligand interaction studies. As the structure of wild type Pu27 is unavailable, modeling its structure as a receptor was a primary goal to understand later the pattern of ligand recognition and ligand binding energetic. The model of the G-quadruplex structure for native Pu27 (5′-TGGGGAGGGTGGGGAGGGTGGGGAAGG-3′), is built over the NMR solution structure of Pu24I¹⁴ using Maestro.³⁷ One thymine and three guanine nucleotides are inserted at the 5′ end and inosine is mutated back to guanine. The cartoon representation of the model built for native sequence is as shown in Fig. 1E, the residue numbers mentioned further are as per the indexing of Fig. 1E.

Modeled-structure is processed with the gas phase energy minimization so as to remove any steric clashes and for the correction of bond orders. The structural integrity of the model is further validated by molecular dynamics simulation studies. Key structural features of Pu24I are well maintained during the simulation run; in addition the 5′ end overhanging bases are capped over the core region of G-tetrads further blocking the probable site for ligand binding i.e. end stacking bases. For studying the binding pattern of ligand we utilized the modeled structure as a starting conformation (and not the end frame of the simulated structure) as the feasible active sites are accessible in it.

2.2. Screening of CHE-like molecules and molecular docking

Based on the similarity search, 14 CHE like molecules are screened from the ‘Supernatural drug compound’ database.³⁸ ADME properties of all the ligands are estimated with the QikProp module of Maestro.³⁹ All molecules have fulfilled the criteria of drug-likeness, their ADME properties are enlisted in Table S1.† Initial binding conformation of CHE and fourteen other CHE like molecules were obtained by molecular docking. 2D structures of ligands were built using maestro and further 3D conformations were generated using Ligprep.⁴⁰ The modeled Pu27 was prepared for molecular docking. GLIDE (Grid based Ligand Docking with energetics) available in Schrödinger software⁴¹ was used for docking. As illustrated by Phan et al., we have considered the 5′ end of G-tetrad core as a ligand binding site.¹⁴ A grid of 10 Å inner-box and 30 Å outer-box was created around the 5′ end of quadruplex structure and the standard precision mode (SP) of docking was used. Docking scores of CHE and 14 CHE like molecules are enlisted in Table S1;† based on the docking score three molecules with better binding were selected for further analysis. They are as shown in Fig. 1, termed as ONE, TWO and THREE. Their docking scores along with the docking score of CHE are enlisted in Table 1.

Table 1 Glide docking score of respective ligands docked on the modeled structure of Pu27

Compound name	Glide docking score^a
a The scores given are approximate values as predicted by the software and not measured experimentally.
CHE	−5.259
ONE	−7.645
TWO	−8.117
THREE	−7.796

---

2.3. Molecular dynamics simulations

Overall 5 systems were considered for simulation analysis i.e. unbound Pu27 and Pu27 in complex with CHE and three CHE like compounds, in the further descriptions these five systems are designated as Pu27 (unbound Pu27), CHE-complex, ONE-complex, TWO-complex and THREE-complex respectively. Same simulation protocol was followed for all the five systems; simulations were performed with the simulation program-AMBER14.⁴² The simple harmonic function used by General Amber Force Field (GAFF) with AM1-BCC charge method was used for parameterization of ligands and ions. Nucleic acid region was parameterized using ff99SB and parmbsc0 force fields of AMBER14.^43,44 Two central K⁺ ions were incorporated into the quartet channels using xleap module of AMBERR14. All systems were neutralized by adding counter-ions and further solvated with 8 Å TIP3PBOX water model.⁴⁵ Standard protocol of AMBER simulation was followed; two step minimization, heating for 50 ps, and equilibration phase for 1 ns. The macroscopic parameters like volume, density, temperature and RMSD of all the systems were equilibrated prior to production run. Equilibrated systems were used for further production run of 50 ns on an NPT ensemble at 300 K temperature and 1 atm pressure, with a step size of 2 fs. Langevin thermostat, barostat were used and SHAKE algorithm was applied.^46–48 Particle mesh Ewald method was applied for long range electrostatic interactions with 0.1 nm grid space of fast Fourier transform grid and the non-bonded cutoff was kept at 12 Å.⁴⁹ Coordinates were recorded in a trajectory at each 10 ps time step. cpptraj module of AMBER, VMD, PyMOL and Chimera were used for the analysis and visualization purposes.^50–53 Simulation data obtained is statistically significant due to long simulation run of 50 ns thus we have considered that the results obtained over the analysis of this data are reliable. Binding free energy of each ligand is estimated over the last 2 ns simulation run of each complex using MMPB/GBSA approach and water density map is generated using cpptraj module. The simulated systems where ligands showed better binding properties were analyzed further for the essential dynamicity in the structure of Pu27 by means of principal component analysis (PCA) and structural properties like Lindemann's coefficient were extracted out from the clustered trajectories. Details of the methodology behind these analyses are described further.

2.4. Binding free energy calculations using MMPB/GBSA method

We followed the standard protocol for calculation of binding free energies. Binding free energy of each ligand was estimated using MMPBSA and MMGBSA approach,^54–58 with the formula:

ΔG_Binding = G_{DNA+ligand complex} − G_DNA − G_ligand

ΔG_Binding = binding free energy, G_{DNA+ligand complex}, G_DNA and G_ligand represent the free energies of respective states.

Free energy of each state was calculated as follows:

G = E_MM + G_PB/GB + G_SA

E_MM = E_νdw + E_ele + E_int

G_SA = γ × SASA + b

E_MM: molecular mechanical energy, G_PB/GB: polar contribution toward solvation energy, Poisson–Boltzmann (PB) or Generalized Born (GB) method used for calculations, G_SA: contributions from nonpolar terms toward solvation energy, E_ele is an electrostatic energy, E_vdw: van der Waals energy, E_int: internal energy (bond, angle and torsional angle energy), SASA: solvent accessible surface area. γ: surface tension proportionality constant (0.0072 kcal mol⁻¹ Å2), b: free energy of nonpolar solvation for a point solute (0 kcal mol⁻¹).

SASA was computed by molsurf using a Linear Combinations of Pairwise Overlaps (LCPO), solute atoms were considered with a probe sphere of 1.4 Å. Binding free energy calculations were averaged over 200 frames taken at the interval of 10 ps over the production run of 2 ns.

2.5. Water-density map analysis method

CPPTRAJ⁵⁰ module of AMBER14 (ref. 42) is used to calculate the occupancy of water molecules around the quadruplex structure for all simulated systems, principle behind this analysis is as follows: occupancy of water is defined as the percentage of time over the trajectory period for which the distance and angle of oxygen atom of water molecule with the corresponding atoms of quadruplex is maintained as present in reference frame. First, the average structure of the 50 ns simulation frames is generated for quadruplex molecule without water molecules and grid of 100 Å with the grid spacing of 0.5 Å is built and fitted to the average structure. Mobility of water is measured over these grid points and intensity of the presence of water molecules at particular grid point gives the water density. For visualization Chimera 1.9 (ref. 53) is used. If particular grid point is populated with water molecules for longer period of time over the entire time scale, that region is covered with dense water contour and conversely if water stays for shorter period the particular region is cover with a lighter water contour which can be examined visually in a XPLOR density map generated in Chimera.

2.6. Principal component analysis and Lindemann's coefficient analysis method

Principal component analysis (PCA) allows the extraction of reference coordinates which elucidate the essential conformational transformations in the macromolecule occurring throughout the simulation run. Thus reduce the vast information obtained by simulation run in to a concise form yet preserve the essence.^59–64 Here we have analyzed the trajectories of CHE-complex, ONE-complex and Pu27 obtained by all atom simulation using PCASuite⁶⁵ package. As the quadruplex structure is showing various phases of conformation during entire simulation run, we divided the trajectories of 50 ns length in to 5 parts each with 10 ns length and further performed the PCA calculations over each of the part individually. Same protocol of PCA is followed for all systems; here procedure for on trajectory is described.

PCA analysis was applied for phosphodiester backbone atoms, trajectory was processed for further analysis using PCAzip module. To remove the oscillations and tumbling motions in the simulated trajectory; all the snapshots were superimposed over the starting frame and after the alignment average coordinates were computed. Further, all the snapshots were superimposed over these average coordinates based on the Gaussian RMSD algorithm⁶⁶ which priorities the weight allotment based on the flexibility of the atom. Covariance matrix of 3N × 3N (N number of atoms) is generated considering coordinates of each atom as a random variable. After the diagonalization of the covariance matrix, relevant eigenvectors were generated which sustained the 90% of the variance. PCZdump module of PCASuite is used to extract out the information from the compressed trajectories such as, eigenvectors with corresponding eigenvalues, hinge points, clustered coordinates of specific eigenvectors and Lindemann's coefficients for each bases. Cumulative percent contribution of eigenvectors in terms of variance was calculated using in house tool. Based on the cumulative percent contribution first two eigenvectors were selected for further analysis as they covered minimum of 20% variance and maximum of 90% variance in some trajectories.

Lindemann's coefficient was calculated for each base for each system, it is taken as an average of 5 sets of calculations for each system. It is calculated for each set of 10 ns trajectory and 5 set of values were averaged and standard errors were calculated for each nucleotide base that too for non backbone atoms. Calculation is based on the formula:^60,67

N = number of atoms, a′ = empirical constant most probable nonbonded near-neighbor distance, ΔR_i² = fluctuation of the atom i.

The porcupine plots were built for visual analysis using modvectors.py code provided by PyMOL and images were built in PyMOL.⁵²

3 Results

3.1. Molecular docking analysis

Analysis of molecular docking gave the superficial idea of possible binding affinities; the docking scores of all the 14 CHE like molecules are enlisted in Table S1.† All the molecules are showing better binding towards Pu27 compared to that of binding of CHE as their docking scores are in more negative range however, we have selected the top three CHE like compounds with better docking scores (ONE, TWO, THREE) for further simulation analysis. Also these molecules have exact core structure of CHE and differ only in the positioning of methoxy groups thus it is likely to explore the influence of position of methoxy group over the binding capacity of respective ligand. Though CHE has a considerable docking score, yet we have channelized it for simulation studies, as experimental data suggests its potential towards the stabilization of Pu27. In addition, it will assist in comparing the binding potential of new molecules with the known potent agent. Docking conformation of selected ligands is described further.

All the four compounds are docked with the feasible conformation as all docking scores are in the negative range (Table 1). CHE and THREE are bound in such a manner that, methoxy groups of ring-A are in polar interactions with the amino group of A15 and aromatic ring formed π–π stacking with G11 (Fig. S1 and S2†). In case of ONE, methoxy groups of ring-A are placed away from the A15 but it is stacked over G11 (Fig. S3†). In case of TWO, ring-D is stacked over G11 and dioxymethyl ring is forming polar interactions with A15 and G13 (Fig. S4†). All the ligands are in plane with the bases A15 and A6, and they are exposed to solvent in major part. Thus we can postulate that like CHE, molecule ONE, TWO and THREE can also bind to Pu27 with an equivalent affinity and similar binding pattern i.e. by means of π–π interactions and polar interactions. Solvent molecules may have role in ligand binding as ligand is extensively exposed to solvent. Further simulation study has been conducted to gain a deeper insight of binding interactions and support our preliminary postulation.

Two out of the three CHE-like molecules (ONE and THREE) which are selected based on the docking score, have been studied as an anti cancer agents. ONE i.e. 12-methoxy chelerythrine is a salt form of plant alkaloid found in Bocconia integrifolia,⁶⁸ THREE i.e. nitidine is a plant alkaloid found in Zanthoxylum nitidium.⁶⁹ Nitidine has been mainly considered for its anti malarial activity,⁷⁰ its anti cancer activity is also well explored.⁷¹ Watanabe et al. have demonstrated that ONE posses anti cancer activity and in MDA-MB-231 cell line it has similar rather slight stronger activity than CHE and THREE.⁷² Bai et al. have performed the detail analysis of sequence selectivity of such plant alkaloids, they have shown that position of methoxy group over ring-A determines the high selectivity of CHE towards GGG repeat sequence over nitidine (THREE) and sanguinarine.⁷³ Thus the interaction patterns of ligands as well as the effect of position of methoxy groups are considered for analysis of simulation data.

3.2. Root mean square displacement analysis

Root mean square displacement (RMSD) with respect to time determines the overall stability of the molecule. Backbone RMSD analysis explains the overall deviation of a structure from its initial conformation. Higher the RMSD more flexible is the structure and if the RMSD is fluctuating it indicates that structure is not stable and yet to achieve the stable most conformation whereas, stable RMSD can be achieved over the simulation run when the conformation converges to a feasible stable state. Here, the ptraj script is used to calculate the RMSD of desired atoms; it performs calculations with the formula:
where, N is a total number of atoms under consideration and δ is the distance between two positions of N pairs of equivalent atoms. As seen in Fig. 2C, all five systems have attained the stable state over the simulation run as all have converged to a steady RMSD values. As elaborated in Fig. 2D, during last 10 ns simulation run all the systems are stabilized around the steady RMSD of ∼4 Å except THREE-complex which is stabilized around ∼6 Å. As the simulation is conducted for long time period of 50 ns the structure of receptor i.e. Pu27 is going through various structural transitions when bound to ligands however in unbound state it has attained a stable conformation within first 10 ns simulation period and sustained the same for rest of the period. Unlike other complexes, TWO-complex also attained stable RMSD values rapidly and it remained in the same state for the rest of the simulation period. In CHE-complex, ONE-complex and THREE-complex overall structure is varying significantly from its initial conformation. The high RMSD fluctuation here should not be misjudged as these ligands are destabilizing the overall structure; visual analyses confirmed that the binding of CHE, ONE and THREE are causing rearrangement in 5′ end bases of Pu27 which is majorly causing the high RMSD in the respective system whereas this phenomenon is completely missing in the dynamics of unbound Pu27 TWO-complex. Thus we can state that TWO has weaker influence over the structure of Pu27 whereas other ligands cause structural alterations. All guanines involved in G-tetrad formation are stable in all the five systems (Fig. 2B). Phan et al. have described the key structural feature of Pu27 such as, presence of G23–A25–G26 triad plane at 3′ end formed by 4^th diagonal loop i.e. ‘G23–A24–A25–G26’ segment,¹⁴ in our simulation this GAG triad is intact in all the five systems (Fig. 2E). RMSD of Loop1 (T10) and Loop3 (T19) are low for all the systems. On the contrary Loop2 formed by residues, G13–G14–A15 are more flexible, Phan et al. also observed the similar phenomenon where G13 and G14 are not well defined in NMR solution spectra may be due to its higher flexibility. Interestingly Loop2 showed higher RMSD in CHE-complex in comparison to that for the others (Fig. 2A). RMSD of ligands (Fig. 2F) is suggesting that all ligands are stable throughout the simulation period with overall RMSD value within the range of ∼0.4 Å to ∼0.6 Å however, THREE shows slight increase in RMSD compared to other ligands which is in correlation to the backbone RMSD data i.e. backbone RMSD of THREE-complex and RMSD of ligand THREE itself are increasing simultaneously. Here simulation is run over long period of 50 ns and all complexes are stable at particular RMSD range for certain time intervals such as, Pu27 is stable for last 40 ns, ONE-complex is stable for last 30 ns, TWO-complex is stable for last 40 ns, CHE-complex is stable for last 40 ns and THREE-complex is at low RMSD for first 20 ns and then stabilized at high RMSD for last 20 ns. In short we can say that binding of THREE is imparting more structural variation in Pu27 as compared to Pu27 alone and its other complex forms however all the structural key features are stable in all the systems.

Fig. 2 RMSD (Å) estimation per unit time (ns), for various fragments of all the five simulated systems; (A) all atom RMSD of Loop2 region (B) all atom RMSD of guanines involved in core stacking (C) backbone RMSD of entire sequence (D) elaboration of backbone RMSD of entire sequence for last 10 ns simulation run (E) all atom RMSD of GAAG capping bases present at the 3′ end region (F) RMSD of ligand molecules.

3.3. B-Factor analysis

RMSD data provides the changes in the overall conformation of the system from its initial state whereas, B-factor offers the fluctuation of each atom or set of atoms over the entire simulation run. It signifies the contribution of each atom in flexibility and stability of the overall structure. B-Factor calculation is performed with cpptraj module of AMBER with the formula:
RMSF represents root mean square fluctuation.

Segment wise B-factor estimation is given in Fig. 3. Similar to RMSD data, all guanines involved in the G-tetrad formation and 3′ end region comprising of GAG triad are inflexible on the other hand 5′ end overhanging bases are highly flexible in all the five systems. Highest flexibility of 5′ bases is observed in THREE-complex, 5′ ends of CHE-complex and ONE-complex possess lower flexible nature than that of THREE-complex, 5′ ends of TWO-complex and Pu27 are least flexible. Loop1 and Loop3 are well defined in the respective NMR solution structure with outward position and Loop2 is ill defined with outward positioning of guanines.¹⁴ In simulation, all loop regions of Pu27 and TWO-complex are lower in flexibility than other complex systems. Loop1 and Loop3 are more flexible in THREE-complex than that of CHE-complex and ONE-complex, Loop2 in CHE-complex is highly flexible in comparison to other complexes. Binding of ligand is majorly influencing the flexibility of active site bases like G5, A15, G16 and G20. Active site bases are not showing much fluctuation in case of Pu27 and TWO-complex. These bases are less flexible in case of CHE-complex and ONE-complex than that of the THREE-complex. Similar to RMSD analyses flexibility patterns also suggest that structures of Pu27 and TWO-complex are alike and binding of TWO has no much influence over the structural pattern of Pu27. Binding of CHE and ONE is inducing dynamicity in the 5′ end overhanging bases and active site residues are getting stabilized. However binding of THREE is somewhere distorting the active site region and arrangement of overhanging bases. Effect of THREE is elaborated in further analysis.

Fig. 3 B-Factor estimation at atomic level, x-axis of each graph represents atom number, y-axis of each graph represents respective B-factor value (A) B-factor of entire sequence of Pu27 in different complexation state, B-factor estimation for different fragments are elaborated further (B) B-factor estimation for Loop1, Loop2 and Loop3 regions (C) B-factor of 5′ end overhanging bases (D) B-factor estimation of bases involved in ligand binding.

3.4. Water density map analysis

Water grid densities have been estimated using cpptraj module of AMBER14, it determines the water structure around the biomolecule. Water occupancy is determined based on the occupancy of oxygen atom of water molecule at particular grid point. Here, the occupancy is defined as the percentage of simulation time at which both the distance and the angle criteria are satisfied. Water density and flexibility of respective region of biomolecules are inversely correlated. As seen in Fig. S5,† CHE-complex and ONE-complex have similar water density, water molecules are structured mainly in the groove regions. In case of TWO-complex and THREE-complex, they are densely solvated and have water structure similar to Pu27. It suggests that binding of CHE and ONE are imparting change in water structure of Pu27 from its unbound state whereas, TWO and THREE do not perturb the dynamics of water in the complex. Number of water molecules around ligands can determine the role of water in stabilizing the binding of ligand with the respective target site. As seen in Fig. S6,† total number of water molecules around the first solvation shell (3.4 Å from ligand) and second solvation shell (5 Å from ligand) of CHE and ONE are increasing as the simulation progressed but in case of TWO and THREE number of water molecules are decreasing. Thus we can state that water molecules are assisting more in the interactions of CHE and ONE with Pu27 than the binding of TWO and THREE. Similar to docked poses all the ligands are highly solvated in the simulated states, these solvent interactions may have role in polar interactions of ligands.

3.5. Conformational behavior of simulated systems

Comparative analyses of the conformational features of Pu27 in unbound state and in complex with the ligand are performed over the simulation trajectories especially, for the key features of Pu24I found in NMR solution structure;¹⁴ they are summarized in Fig. 4. G23–A25–G26 triad maintains the stable planar arrangement in all the five systems and interaction pattern is same as in the original structure of Pu24I i.e. G23–G26 are paired with the Hoogsteen hydrogen bonding, G23–A25 are interacting with GA sheared hydrogen bonding and, A24 is showing π–π stacking with G23. The out ward conformations of Loop1 (T10 in Pu27) and Loop3 (T19 in Pu27) as found in Pu24I, are conserved in all the five systems whereas, conformation of Loop2 (G13, G14, A15 in Pu27) is different in all the systems. In Pu24I, G13 and G14 are positioned outward with no stacking interaction and A15 is stabilized by mismatch interaction with A6 along with stacking over the first G-plane. In both Pu27 and TWO-complex, G13 and G14 are facing out and form π–π stacking interactions among themselves, they also stack with the T1 i.e. first residue of 5′ end thus facilitating the capping of 5′ end overhang over the G-tetrad core. In both CHE-complex and ONE-complex, G13 and G14 are facing out with stacking interactions among themselves, however no stacking interactions with T1 was found and capping of 5′ end overhang over the G tetrad is missing. In case of THREE-complex, G13 and G14 are facing inward with stacking among themselves, stacking interaction with T1 and capping of 5′ end overhang over the G tetrad is also lacking. The planar conformation of A6–A15 stacked over the first G-plane is stable in all the systems, in CHE-complex and ONE-complex these bases are offering π–π stacking interactions for ligand binding.

Fig. 4 Conformational behavior of Pu27 in different binding states: GAG triad at 3′ end found to be stable in all the five systems also; thymine bases in Loop1 and Loop3 are facing outward in all the five systems. Loop2 bases and ligand binding site near A–A mismatch region are showing conformational variation with change in complexation state.

The conformational arrangement described earlier for the 5′ end overhanging bases is based on the average geometry of the respective structure however; these are attained through the dynamic transformations of the starting conformation. In unbound state i.e. Pu27, the outward stacked conformations of Loop2 and capping of 5′ end overhanging bases are attained rapidly but binding of ligand delays this rearrangement. In presence of ONE, when ligand is localized out of the binding site the 5′ end overhanging bases and Loop2 slowly attained the conformation similar to that of the unbound state. The binding pattern and re-localization of ONE from outer side to inner pocket is facilitated by means of characteristic crawling of methyl fragments of methoxy groups over the aromatic rings of nucleotide bases via CH₃–π interactions. Methoxy group at the 8th position of ring-A initiated the CH₃–π interactions with G11 and crawled inside leading to CH₃–π stacking interactions of methoxy group at 7^th position with G11. Methoxy group at the 8^th position further perturbed the initial pose of G5 (guanine present in the 5′ end overhang) thus destructing the stacking between T1 and G14. It further instigated the CH₃–π stacking with G5 and crawled inside the space between A6 and G5, eventually ONE is intercalated between the A6–A15 plane and G5. CHE interaction pattern is differing slightly than that of the binding of ONE. Binding of CHE is rapid, even before the assembly of 5′ end overhanging bases CHE is intercalated between A6–A15 plane and G5 further restricting the movement of overhanging bases thus capping is adjourned. In the final conformation of CHE, methoxy groups at 7^th and 8^th position of ring-A are localized in a way that they captured the backbone region of G5 and A6 thus arresting their movement. In the final binding mode of both CHE and ONE, at both the faces of ring-A and ring-B they are forming π–π stacking interactions with the aromatic rings of A6 and G5, ring-C and ring-D are stacked over the aromatic ring of A15 further more dioxymethyl ring is in interaction with the backbone atoms of A15. Unlike the binding pattern of CHE and ONE, TWO is localized outside the binding site with the deficiency in CH₃–π stacking interactions and other polar interactions. However, the 5′ end capping rearrangement is retained in TWO-complex but the process is delayed as compared to the unbound state. Binding pattern of THREE is found to be indistinct, for major simulation span THREE is out of the binding site with no appropriate stacking however of 5′ end bases are not capped over the core region of G-tetrads. Towards the end of the simulation run, 5′ end overhanging bases are distorted and THREE is intercalated between G5 and G4. Movies of simulated trajectories are provided in the ESI† for visual confirmations. Uncertainty in the interaction patterns of TWO and THREE is also reflected in the estimation of binding energy with the MMPBSA approach; ΔG_Binding estimated for TWO and THREE have higher standard deviations and standard errors i.e. binding energy is fluctuating over the simulation run hence, they are not considered here for analyzing the thermodynamic parameters.

3.6. Binding free energy estimation

As per the performance based analysis of MMPBSA and MMGBSA methods it has been stated that; binding free energy prediction is closer to absolute values with MMPBSA calculations however, MMGBSA results show better correlation with the experimental findings.⁷⁴ Thus ranking of the binding affinities among different ligands is possible with the MMGBSA approach whereas absolute binding energy estimation is feasible with the MMPBSA approach. Here, we have discussed the binding free energy estimation of CHE and ONE which are calculated with both the methods, details of energetic contributions are as enlisted in Table 2. Binding of these ligands is preferable as the overall binding free energies are in negative range with the equivalent values (MMPBSA method: ∼−18 kcal mol⁻¹). As per MMGBSA estimations, ONE possesses slight higher binding affinity than that of CHE and these results are in correlation with the inhibitory properties observed in MDA-MB-231 cell line⁷² (discussed earlier). Electrostatic interactions are showing major contribution in binding energies however they are counterbalanced by polar solvation energies. Ultimately, van der Waals interactions are determining the binding affinity. Here, additional methoxy group at 12^th position of ONE is boosting the binding affinity over that of CHE by means of enhanced van der Waals interactions.

Table 2 MMPB(GB)SA calculations of binding free energy components of ligand ONE and TWO when bound to Pu27, analyzed over last 2 ns simulation period^a

Energy^b (kcal mol^−1)	CHE	ONE
a EVDWAALS: non-bonded van der Waals energy, EEEL: non-bonded electrostatic energy, EPB: polar solvation energy (PB), EGB: polar solvation energy (GB), ENPOLAR: non polar solvation energy from repulsive solute–solvent interactions (PB), ESURF: non polar solvation energy (GB), EDISPER: non-polar contribution to solvation energy from attractive solute–solvent interactions (PB), ΔGGas: EVDWAALS + EEEL + internal energy, ΔGSolv: polar solvation energy + non polar solvation energy, ΔGBinding energy: binding free energy (ΔGGas + ΔGSolv).b The energies given are approximate values as predicted by the software and not measured experimentally.
EVDWAALS	−41.54 ± 2.19	−46.16 ± 2.25
EEEL	−606.25 ± 9.14	−600.28 ± 13.93
EPB	611.64 ± 9.44	607.77 ± 14.62
EGB	620.59 ± 9.55	615.33 ± 14.84
ENPOLAR	−21.33 ± 0.95	−22.46 ± 1.21
ESURF	−2.37 ± 0.12	−2.40 ± 0.21
EDISPER	38.72 ± 0.98	42.95 ± 1.57
ΔGGas	−647.79 ± 10.08	−646.44 ± 15.12
ΔGSolvPB	629.03 ± 9.39	628.26 ± 15.11
ΔGSolvGB	618.22 ± 9.53	612.93 ± 14.68
ΔGBinding energyPB	−18.76 ± 2.74	−18.18 ± 3.15
ΔGBinding energyGB	−29.57 ± 2.15	−33.51 ± 2.17

---

3.7. Principal component analysis

Principal component analysis (PCA) also termed as essential dynamics has been applied in the analysis of conformational behavior of biomolecules.^59–64 Haider et al. have determined the dynamic behavior of multimers of telomeric G-quadruplex using PCA analysis, they have found that loop regions are highly flexible however stacked guanines are rigid and ligand binding enhance the rigidity of the quadruplex.⁵⁹ In our study, CHE and ONE are found to be effective stabilizing agents of Pu27 thus we further explored the effect of CHE and ONE over the dynamicity of Pu27 in comparison to that of its unbound form. The binding of CHE and ONE have induced the changes in the overall conformation of Pu27 which have been reflected in the high RMSD values of these simulated systems. To remove the noise and to determine the major dynamic regions behind these alterations we have performed the PCA of MD trajectories of unbound Pu27, CHE-complex and ONE-complex using PCASuite package.⁶⁵

As seen in Fig. 5, first two eigenvectors of CHE-complex and ONE-complex and in case of unbound Pu27 first 6 eigenvectors cover >50% of variance. Here we have analyzed the dynamicity covered by first two eigenvectors of these systems. Projections of the porcupine plot denotes the dynamicity of the particular backbone phosphorous atom however direction is not considered as the visual analysis of compressed trajectories suggested the motion of flexible regions is back and forth along the axis of the respective projection. The porcupine plot of unbound Pu27 for first two eigenvectors (Fig. 6A and F) suggest that for first 10 ns simulation run the 5′ end over hanging bases are highly flexible and their dynamicity aids in the aromatic interaction of T1 with G14 thus forming cap over the 5′ end, on the other hand other bases are rigid. For next 30 ns (Fig. 6B–D and G–I), the dynamicity of 5′ end overhanging bases is reduced as T1, G14, G13 are stacked over each other. During last 10 ns simulation run the bases of Loop2 (G13, G14) and T1 are moving all together without breaking the aromatic interactions.

Fig. 5 ‘Cumulative% contribution of each eigenvector’ for each trajectory derived from PCA analysis of (A) unbound Pu27 system, (B) CHE-complex and (C) ONE-complex. 1: 0–10 ns, 2: 11–20 ns, 3: 21–30 ns, 4: 31–40 ns, 5: 41–50 ns time frames of each system are plotted.

Fig. 6 Porcupine plots of first (A–E) and second eigenvectors (F–J) of unbound Pu27. Time frames are mentioned for respective plots (images are prepared in PyMOL).

If compared the effect of binding of CHE and ONE over the dynamic behavior of Pu27 (Fig. 7 and S7†) with that of the unbound state (Fig. 6), it is clear that binding of ligand enhance the flexibility of the 5′ overhanging bases. In case of CHE-complex (Fig. 7), CHE is entering in to the intercalation site during first 20 ns time steps and locking the overhanging bases away from Loop2 thus preventing the cap formation. In case of ONE-complex ligand is outside the pocket for first 10 ns simulation run thus like unbound state of Pu27, 5′ overhanging bases are involved in the cap formation with Loop2 (Fig. S7A and F†). During 11–20 ns, the motion of ONE is coordinated with the movement of G5 (Fig. S7B and G†) thus facilitating the intercalation of ONE. For last 20 ns the conformation of 5′ overhanging bases is same as that of CHE-complex i.e. away from Loop2. In all the three systems guanines involved in core stacking regions are rigid.

Fig. 7 Porcupine plots of first (A–E) and second eigenvectors (F–J) of CHE-complex. Time frames are mentioned for respective plots. Cyan color is for quadruplex structure and red color is for chelerythrine molecule (images are prepared in PyMOL).

3.8. Lindemann's coefficient analysis

Zhou et al. have described the application of Lindemann's criterion in the molecular dynamics simulation of protein molecules, they have determined that residues which are buried inside the core of protein are more solid-like while the residues located over the surface area are more liquid-like.⁶⁷ The buried residues with low Lindemann's coefficient (<0.15) have solid nature (less flexible) and provide the stability while surface residues with high Lindemann's coefficient (>0.15) have liquid nature provide the conformational flexibility needed for activity of respective protein, Jamroz et al. have determined the flexibility of proteins based on the Lindemann's coefficient.⁶⁰ Here we have utilized the same concept to determine the solid–liquid nature of the nucleotide bases using PCASuite package,⁶⁵ calculation details are provided in the ESI.†

As seen in Fig. 8, the guanine bases involved in core stacking (residue number: 7–9, 11, 12, 16–18, 20–22, 27) are in solid state in all the three systems nonetheless binding of ONE and CHE further reduce the Lindemann's coefficient for the same bases in comparison to that of the unbound state it indicates that binding of the ligand further stabilizes these guanines. Thymine bases in Loop1 (T10) and Loop3 (T19) are in liquid state in all the three systems however binding of ligands slightly reduce the flexibility of the same in comparison to the unbound state. Loop2 bases (G13, 14) are in liquid state, G13 has high standard error thus we can state that flexibility of this base is unaltered upon ligand binding whereas G14 is found to be highly flexible in CHE-complex than other two systems. A6 and A15 involved in mismatch pairing are equally stable in all the systems. The stretch of G23–A–A–G26 nucleotides which protect 3′ end are more solidified upon ligand binding; binding of ONE has slightly more stabilizing effect than that of the CHE. Here we find that flexibility of 5′ end overhanging bases (residue number: 1–5) is enhanced over the binding of CHE and ONE.

Fig. 8 Lindemann's coefficient per residues is determined for clustered trajectories of unbound Pu27, CHE-complex and ONE-complex. Coefficients are calculated for all the five trajectories separately and average is plotted along with the standard error bars. Value < 0.15 indicates solid nature of residue whereas value > 0.15 indicates liquid nature.

4 Discussion and conclusion

The structure of the entire sequence of Pu27 is tedious to obtain, as in solution state it interchanges between the two topological arrangements though only one conformation is found to be involved in biological activities. Thus as far now, structural analysis of ligand induced stabilization of Pu27 is carried out either with the altered sequence or with the truncated form of native sequence. Ghosh et al. have offered a novel mechanism of action of CHE for its anti cancer activity i.e. ‘CHE stabilizes the secondary structure of G-quadruplex formed in the promoter region of c-MYC oncogene (Pu27) and further arrests it's expression’. Thus here we have performed the interaction studies of CHE and three CHE like natural compounds with the unaltered native sequence of Pu27, aiming to (i) explore the atomistic details of binding pattern of CHE and conformational behavior of Pu27 upon binding of ligand, (ii) unravel the influence of overhanging 5′ end bases over the binding of ligands to Pu27 and (iii) determine the position specific role of methoxy groups in the CHE skeleton modulating the binding ability of CHE and CHE like molecules.

All the key structural features observed in previous reported structures of Pu27 formed of either altered or truncated sequences are conserved in the simulation of the modeled structure of native sequence thus it can be utilized for the study purposes. In simulation it is observed that, the overhanging bases at 5′ end of unbound state are folding over the core region of G-tetrad whereas, in presence of ligands this capping rearrangement is either restricted or delayed. Binding of CHE and ONE restricts the capping of 5′ end overhanging bases instead they gain extended conformation, binding of TWO delays this rearrangement whereas binding of THREE completely distorts the 5′ end. The CH₃–π interactions are playing key role in the binding phenomenon of ligands; presence of methoxy groups at 7^th and 8^th position of ring-A is found to be crucial for the characteristic crawling motion of ligand to dive into the binding site in Pu27. In case of CHE and ONE rightly positioned methoxy groups assisted them to ultimately intercalate between G5 (part of overhanging bases) and A6–A15 plane. Methoxy group at 12^th position of ONE is further enhancing the binding affinity towards Pu27 by raising the energy contribution through van der Waals interactions. 12^th position of the benzo[c]phenanthridinium compounds has been found to play crucial role in determining their cytotoxic properties.^72,75–78 Mackay et al. have reported the influence of 12-alkoxy group of N-methyl benzo[c]phenanthridinium compounds over their anti cancer activity. They found that non functionalized shorter aliphatic chains as 12-alkoxy group show better anticancer activity.⁷⁶ Lynch et al. have demonstrated that 12-ethoxy substitution over N-methyl benzo[c]phenanthridinium compounds enhances their anti cancer activity as well as it imparts the sequence specificity.⁷⁷ It has been observed that incorporation of methoxy group at 12^th position enhances the anticancer activity of macarpine⁷⁸ also, improves its DNA binding.⁷⁵ Similarly our observation suggests that substitution of methoxy group at 12^th position of CHE is enhancing the binding affinity towards Pu27.

Intercalated ligands are stabilized majorly through π–π stacking interactions with aromatic rings of respective bases. The ligand mediated restricted motion of G5 and its localization over the core region of G-tetrads prevents the rearrangement of 5′ end overhanging bases and thus capping. In case of TWO, methoxy groups are placed at 9^th and 10^th position of ring-A; this shift hinders the desired CH₃–π interactions and restricts the entry of ligand in to the binding site thus the overall conformation of G-quadruplex remains unperturbed from that of the unbound state. In case of THREE these methoxy groups are located over 8^th and 9^th position of ring-A, which hampers the CH₃–π interactions and prevents its proper binding. Ghosh et al. have confirmed that two molecules of CHE binds to Pu27 by means of end-stacking;¹⁵ the 3′ end is enclosed with the GAG triad and groove regions are not accessible due to high flexibility of loop fragments thus, the plausible binding site will be the 5′ end of G-tetrad core. In simulation we found that, in unbound state the ligand binding site at 5′ end is concealed by the cap formed by the overhanging bases however, binding of CHE unwraps this site as it intercalates within the end stacking bases (A6–A15, G5). CHE further arrests the dynamicity of overhanging bases thus the binding site is retained in an open state. This may be facilitating the binding of second molecule of CHE over the 5′ end i.e. stacking over G7–G11/G16–G20 region. Such dual mode of binding is also feasible for ONE as it shares the similar binding characteristics of CHE.

In summary, CHE and ONE (12-methoxy CHE) intercalate within the plane of A6–A15 and G5 at the 5′ end of G-quadruplex (secondary structure of Pu27) majorly through CH₃–π and π–π stacking interactions and may assist in the binding of second molecule. G5 assists the binding of ligands, methoxy groups at 7^th and 8^th position of ring-A are essential for proper binding, methoxy group at 12^th position enhances the binding affinity. These key findings can act as rational behind the designing of selective leads for stabilization of Pu27.

Conflict of interest

Authors declare no conflict of interests.

Abbreviations

CHE	Chelerythrine
Pu27	Unbound Pu27
ONE	12-Methoxy chelerythrine (database ID: 4386-SN00288441)
TWO	(Database ID: 4386-SN00375322)
THREE	Nitidine (database ID: 4386-SN00155937)

Acknowledgements

SC thanks Department of Science and technology for Ramanujan Fellowship (SR/S2/RNJ-78/2010 (G)). JB thanks the DBT fellowship from the project (BT/PR6627/GBD/27/440/2012).

References

1. C. A. Spencer and M. Groudine, Adv. Cancer Res., 1991, 56, 1–48.
2. D. M. Miller, S. D. Thomas, A. Islam, D. Muench and K. Sedoris, Clin. Cancer Res., 2012, 18, 5546–5553.
3. S. Pelengaris, B. Rudolph and T. Littlewood, Curr. Opin. Genet. Dev., 2000, 10, 100–105.
4. K. B. Marcu, S. A. Bossone and A. J. Patel, Annu. Rev. Biochem., 1992, 61, 809–858.
5. L. M. Facchini and L. Z. Penn, FASEB J., 1998, 12, 633–651.
6. M. Cooney, G. Czernuszewicz, E. H. Postel, S. J. Flint and M. E. Hogan, Science, 1988, 241, 456.
7. T. Simonsson, M. Pribylova and M. Vorlickova, Biochem. Biophys. Res. Commun., 2000, 278, 158–166.
8. T. Simonsson, M. Kubista and P. Pecinka, Nucleic Acids Res., 1998, 26, 1167–1172.
9. T. Tomonaga and D. Levens, Proc. Natl. Acad. Sci. U. S. A., 1996, 93, 5830–5835.
10. E. F. Michelotti, T. Tomonaga, H. Krutzsch and D. Levens, J. Biol. Chem., 1995, 270, 9494–9499.
11. E. H. Postel, S. J. Berberich, J. W. Rooney and D. M. Kaetzel, J. Bioenerg. Biomembr., 2000, 32, 277–284.
12. A. Siddiqui-Jain, C. L. Grand, D. J. Bearss and L. H. Hurley, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 11593–11598.
13. H.-J. Kang and H.-J. Park, J. Comput.-Aided Mol. Des., 2015, 29, 339–348.
14. A. T. n. Phan, V. Kuryavyi, H. Y. Gaw and D. J. Patel, Nat. Chem. Biol., 2005, 1, 167–173.
15. S. Ghosh and D. Dasgupta, Biochem. Biophys. Res. Commun., 2015, 459, 75–80.
16. T.-M. Ou, Y.-J. Lu, C. Zhang, Z.-S. Huang, X.-D. Wang, J.-H. Tan, Y. Chen, D.-L. Ma, K.-Y. Wong and J. C.-O. Tang, J. Med. Chem., 2007, 50, 1465–1474.
17. R. Rocca, F. Moraca, G. Costa, S. Alcaro, S. Distinto, E. Maccioni, F. Ortuso, A. Artese and L. Parrotta, Molecules, 2014, 20, 206–223.
18. X. Ji, H. Sun, H. Zhou, J. Xiang, Y. Tang and C. Zhao, Nucleic Acid Ther., 2012, 22, 127–136.
19. S. A. Ohnmacht and S. Neidle, Bioorg. Med. Chem. Lett., 2014, 24, 2602–2612.
20. J. Jana, R. K. Kar, A. Ghosh, A. Biswas, S. Ghosh, A. Bhunia and S. Chatterjee, Mol. BioSyst., 2013, 9, 1833–1836.
21. R. K. Kar, P. Suryadevara, J. Jana, A. Bhunia and S. Chatterjee, J. Biomol. Struct. Dyn., 2013, 31, 1497–1518.
22. S. Mondal, J. Bhat, J. Jana, M. Mukherjee and S. Chatterjee, Mol. BioSyst., 2016, 12, 18–22.
23. I. Bessi, C. Bazzicalupi, C. Richter, H. R. A. Jonker, K. Saxena, C. Sissi, M. Chioccioli, S. Bianco, A. R. Bilia and H. Schwalbe, ACS Chem. Biol., 2012, 7, 1109–1119.
24. S. Rybáková, M. Rájecký, J. Urbanová, K. Pěnčíková, E. Táborská, R. Gargallo and P. Táborský, Chem. Pap., 2013, 67, 568–572.
25. S. Ghosh, A. Kar, S. Chowdhury and D. Dasgupta, Biochemistry, 2013, 52, 4127–4137.
26. Y. Ma, T.-M. Ou, J.-Q. Hou, Y.-J. Lu, J.-H. Tan, L.-Q. Gu and Z.-S. Huang, Bioorg. Med. Chem., 2008, 16, 7582–7591.
27. K. Bhadra and G. S. Kumar, Biochim. Biophys. Acta, 2011, 1810, 485–496.
28. M. Franceschin, L. Rossetti, A. D'Ambrosio, S. Schirripa, A. Bianco, G. Ortaggi, M. Savino, C. Schultes and S. Neidle, Bioorg. Med. Chem. Lett., 2006, 16, 1707–1711.
29. S. Ghosh, S. K. Pradhan, A. Kar, S. Chowdhury and D. Dasgupta, Biochim. Biophys. Acta, 2013, 1830, 4189–4201.
30. W.-J. Zhang, T.-M. Ou, Y.-J. Lu, Y.-Y. Huang, W.-B. Wu, Z.-S. Huang, J.-L. Zhou, K.-Y. Wong and L.-Q. Gu, Bioorg. Med. Chem., 2007, 15, 5493–5501.
31. S. K. Pradhan, D. Dasgupta and G. Basu, Biochem. Biophys. Res. Commun., 2011, 404, 139–142.
32. J. M. Herbert, J. M. Augereau, J. Gleye and J. P. Maffrand, Biochem. Biophys. Res. Commun., 1990, 172, 993–999.
33. S. Ghosh, J. Jana, R. K. Kar, S. Chatterjee and D. Dasgupta, Biochemistry, 2015, 54, 974–986.
34. A. T. n. Phan, Y. S. Modi and D. J. Patel, J. Am. Chem. Soc., 2004, 126, 8710–8716.
35. S. Burge, G. N. Parkinson, P. Hazel, A. K. Todd and S. Neidle, Nucleic Acids Res., 2006, 34, 5402–5415.
36. R. I. Mathad, E. Hatzakis, J. Dai and D. Yang, Nucleic Acids Res., 2011, 39, 9023–9033.
37. Schrödinger Release 2013-3: Maestro, version 9.6, Schrödinger, LLC, New York, NY, 2013.
38. M. Dunkel, M. Fullbeck, S. Neumann and R. Preissner, Nucleic Acids Res., 2006, 34, D678–D683.
39. Small-Molecule Drug Discovery Suite 2013-3: QikProp, version 3.8, Schrödinger, LLC, New York, NY, 2013.
40. Schrödinger Release 2013-3: LigPrep, version 2.8, Schrödinger, LLC, New York, NY, 2013.
41. Small-Molecule Drug Discovery Suite 2013-3: Glide, version 6.1, Schrödinger, LLC, New York, NY, 2013.
42. D. A. Case, V. Babin, J. Berryman, R. M. Betz, Q. Cai, D. S. Cerutti, T. E. Cheatham Iii, T. A. Darden, R. E. Duke and H. Gohlke, AMBER-14.
43. J. Wang, W. Wang, P. A. Kollman and D. A. Case, J. Mol. Graphics Modell., 2006, 25, 247–260.
44. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
45. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
46. R. J. Loncharich, B. R. Brooks and R. W. Pastor, Biopolymers, 1992, 32, 523–535.
47. J. A. Izaguirre, D. P. Catarello, J. M. Wozniak and R. D. Skeel, J. Chem. Phys., 2001, 114, 2090–2098.
48. V. Kräutler, W. F. van Gunsteren and P. H. Hünenberger, J. Comput. Chem., 2001, 22, 501–508.
49. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
50. D. R. Roe and T. E. Cheatham Iii, J. Chem. Theory Comput., 2013, 9, 3084–3095.
51. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
52. W. L. DeLano, The PyMOL molecular graphics system, 2002.
53. E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, D. M. Greenblatt, E. C. Meng and T. E. Ferrin, J. Comput. Chem., 2004, 25, 1605–1612.
54. R. N. Shinde and M. E. Sobhia, J. Mol. Graphics Modell., 2013, 45, 98–110.
55. T. E. Cheatham Iii, J. Srinivasan, D. A. Case and P. A. Kollman, J. Biomol. Struct. Dyn., 1998, 16, 265–280.
56. B. Kuhn and P. A. Kollman, J. Med. Chem., 2000, 43, 3786–3791.
57. B. Honig and A. Nicholls, Science, 1995, 268, 1144–1149.
58. A. Onufriev, D. Bashford and D. A. Case, Proteins: Struct., Funct., Bioinf., 2004, 55, 383–394.
59. S. Haider, G. N. Parkinson and S. Neidle, Biophys. J., 2008, 95, 296–311.
60. M. Jamroz, M. Orozco, A. Kolinski and S. Kmiecik, J. Chem. Theory Comput., 2012, 9, 119–125.
61. A. Amadei, A. Linssen and H. J. C. Berendsen, Proteins: Struct., Funct., Bioinf., 1993, 17, 412–425.
62. A. Amadei, A. B. M. Linssen, B. L. De Groot, D. M. F. Van Aalten and H. J. C. Berendsen, J. Biomol. Struct. Dyn., 1996, 13, 615–625.
63. D. M. F. Van Aalten, J. B. C. Findlay, A. Amadei and H. J. C. Berendsen, Protein Eng., 1995, 8, 1129–1135.
64. H. Yamaguchi, D. M. F. Van Aalten, M. Pinak, A. Furukawa and R. Osman, Nucleic Acids Res., 1998, 26, 1939–1946.
65. T. Meyer, C. Ferrer-Costa, A. Perez, M. Rueda, A. Bidon-Chanal, F. J. Luque, C. A. Laughton and M. Orozco, J. Chem. Theory Comput., 2006, 2, 251–258.
66. K. L. Damm and H. A. Carlson, Biophys. J., 2006, 90, 4558–4573.
67. Y. Zhou, D. Vitkup and M. Karplus, J. Mol. Biol., 1999, 285, 1371–1375.
68. S. M. Oechslin, G. M. König, K. Oechslin-Merkel, A. D. Wright, A. D. Kinghorn, O. Sticher and M. Miyagawa, J. Nat. Prod., 1991, 54, 519–524.
69. H. R. Arthur, W. H. Hui and Y. L. Ng, J. Chem. Soc., 1959, 1840–1845.
70. J. Bouquet, M. Rivaud, S. Chevalley, E. Deharo, V. Jullian and A. Valentin, Malar. J., 2012, 11, 67.
71. M. Del Poeta, S.-F. Chen, D. Von Hoff, C. C. Dykstra, M. C. Wani, G. Manikumar, J. Heitman, M. E. Wall and J. R. Perfect, Antimicrob. Agents Chemother., 1999, 43, 2862–2868.
72. T. Watanabe, Y. Ohashi, R. Yoshino, N. Komano, M. Eguchi, S. Maruyama and T. Ishikawa, Org. Biomol. Chem., 2003, 1, 3024–3032.
73. L.-P. Bai, Z.-Z. Zhao, Z. Cai and Z.-H. Jiang, Bioorg. Med. Chem., 2006, 14, 5439–5445.
74. T. Hou, J. Wang, Y. Li and W. Wang, J. Comput. Chem., 2011, 32, 866–877.
75. M. Rájecký, K. Šebrlová, F. Mravec and P. Táborský, PLoS One, 2015, 10, e0129925.
76. S. P. Mackay, L. Comoe, B. Desoize, O. Duval, J.-C. Jardillier and R. D. Waigh, Anti-Cancer Drug Des., 1998, 13, 797–813.
77. M. A. Lynch, O. Duval, A. Sukhanova, J. Devy, S. P. MacKay, R. D. Waigh and I. Nabiev, Bioorg. Med. Chem. Lett., 2001, 11, 2643–2646.
78. T. Ishikawa, T. Saito and H. Ishii, Tetrahedron, 1995, 51, 8447–8458.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra04671a

---