Computational DNA binding studies of (–)-epigallocatechin-3-gallate

Abstract

The catechin family of molecules that are present in the leaves of green tea has been under investigation since the antioxidant and anti-inflammatory properties of tea were discovered. Among multiple proposed therapeutic targets of these molecules, the direct interaction with nucleic acids has been proposed and experimentally observed but without clear knowledge about the potential binding modes between these ligands and DNA. One of these catechin structures, (–)-epigallocatechin gallate (EGCG), has three aromatic rings that could interact with double-stranded DNA via terminal base-pair stacking, intercalation, or through groove binding. Using enhanced sampling techniques and molecular dynamics simulations, we have found a stable complex between the EGCG ligand and DNA through intercalation of the trihydroxybenzoate aromatic ring and an ApC step. Moreover, we have calculated the absorption spectra of four possible binding modes and compared these to absorption profiles reported in the literature, and explored the possible DNA sequence preference for the EGCG ligand to bind. Our results suggest that an intercalative mode of interaction through the major groove is possible between the EGCG ligands and DNA with apparently very little DNA sequence selectivity.

Keywords:

• epigallocatechin gallate
• DNA interactions
• AMBER
• molecular dynamics

Introduction

Antioxidant molecules aid in the reduction of chronic disease and are commonly found in fruits and vegetables (Weisburger). Molecules that contain three or more hydroxyl groups in aromatic rings, known as polyphenols, have been shown to possess antioxidative properties. A sub-group of polyphenols consists of a series of catechin molecules that include (–)-epigallocatechin gallate (EGCG), (–)-epigallocatechin (EGC), (–)-epicatechin gallate (ECG), and (–)-epicatechin (EC). These compounds are found in the leaves of Camellia sinensis (tea), and as early as the 1980s their anti-inflammatory and antioxidant properties were evident (Chen et al., 2011).Anti-cancer properties of these polyphenols, mainly of EGCG (Du et al., 2012), were first evidenced by Want and collaborators using in vivo studies in an animal model (Wang et al., 1995). This opened a door to perform further testing and therapeutic assessment of catechin-type molecules found in tea extracts. The effectiveness of EGCG as a tumor suppressor has been routinely observed in hundreds of studies (Bansal, Syan, Mathur, & Choudhary, 2011; Beltz, Bayer, Moss, & Simet, 2006; Chen et al., 2011; Kazi et al., 2002; Rahmani, Al Shabrmi, Allemailem, Aly, & Khan, 2015; Shankar, Ganapathy, & Srivastava, 2007). Progress has been made to understand the mechanism of action for EGCG, which can lead to a limited amount of tumor regression that involves multiple sites and multiple cellular pathways, including the inhibition of kinases, angiogenesis inhibition, antioxidant activity, and modulation of transcription factors (Rahmani et al., 2015). The interaction with nucleic acids is then a main focus of research to understand the effectiveness and mechanism of action of EGCG. This is supported by mass spectroscopy studies that have shown direct interaction of EGCG with double-stranded DNA and telomeric quadruplex DNA (Mikutis et al., 2013).

Despite the many studies and publications, the reality of a positive pharmacological effect of the EGCG molecule is unsure. In vitro studies using cancer cell lines suggest that natural antioxidant supplements may actually aid in promoting the growth of cancerous cells rather than to prevent it (Lu, Ou, & Lu, 2013). Regardless, studies show that EGCG inhibits cell growth and aids in promoting apoptosis, and known potent anti-cancer drugs that interact directly with DNA have been observed to enhance their effect when EGCG is included. The compound cisplatin, which is a transition metal drug routinely included in chemotherapy, showed increased activity when EGCG was included (Chan, Soprano, Weinstein, & Fong, 2006; Wang et al., 2015), and restored the toxicity to previously desensitized lung cancer cells (Zhang, Wang, Han, Zhou, & Sun, 2015). Similar synergistic effects were observed with oxaliplatin, a cisplatin derivative (Hu et al., 2015). Cooperative effect was observed when the anti-cancer, DNA binding compound bleomycin was used in conjunction with EGCG, showing a potentiated anti-proliferative effect in pancreatic cancer cells (Bimonte et al., 2015). Additionally, interaction of EGCG in single-strand DNA chains showed increased nuclease activity in whole blood and isolated lymphocytes cells (Yuquan, Takeshita, & Morimoto, 2001). Using surface plasmon resonance assays and cold spray ionization-mass spectrometry, Fujiki and collaborators found that EGCG binds to DNA and RNA polymers between 1:1 and 3:1 stoichiometry (Kuzuhara, Sei, Yamaguchi, Suganuma, & Fujiki, 2006). They additionally detected that multiple binding conformations are present in EGCG–DNA and EGCG–RNA complexes. Electrochemical studies of EGCG when in contact with fish testes DNA revealed a direct interaction with DNA via possible intercalation between the DNA base-pairs (Zheng, Chen, Hoshi, Anzai, & Li, 2006). Dasgupta and collaborators performed the analysis on the interaction of EGCG in the presence of copper ions and calf-thymus DNA using several spectroscopic techniques (Ghosh, Sahoo, Jana, & Dasgupta, 2008). They observed the fluorescence quenching of ethidium bromide bound to DNA with increasing concentrations of EGCG, which implies that EGCG molecules may interact with double-stranded DNA through the insertion of one of its rings between base-pairs.

Spectroscopic assays that are routinely used to study the interaction of small ligands to DNA include UV-visible absorption and fluorescence quenching (Sirajuddin, Ali, & Badshah, 2013). In UV-vis experiments, the changes in the position of the maximum value of the absorption provide information on whether DNA-ligand interactions are occurring. The literature suggests that a hypochromic effect with no absorption wavelength shift is related to an intercalative mode of interaction, whereas a hyperchromic effect combined with shifting in the absorption wavelength (either blue or redshift) points to a groove binding type of interaction (Sirajuddin et al., 2013). In lack of a solution structure of the actual EGCG + DNA complex, spectroscopic information currently only provides evidence that an actual interaction between ligand and DNA is occurring, but the actual mode of interaction is speculative at best, since multiple experiments found in the literature point to contradictory experimental observations.

In order to study the modes of interaction and stability of EGCG bound to DNA, we performed molecular dynamics (MD) simulations in which each of the aromatic rings of the ligand structure (A, B and D, Figure 1) is manually intercalated between an AT base-step. The intercalation process was studied using equilibrium umbrella sampling (US) with the aid of a biasing force to drive the system between the desired states. The use of this enhanced sampling technique has been used before to study the intercalation process of several DNA binders (Galindo-Murillo, Garcia-Ramos, Ruiz-Azuara, Cheatham, & Cortes-Guzman, 2015) and is a convenient methodology to study a particular problem that otherwise would require millisecond sampling rate, which is still not reachable without the use of specialized MD-specific hardware, such as the ANTON supercomputer (Shaw et al., 2009).

Figure 1. Structure of (–)-epigallocatechin-3-gallate (EGCG). Rings A and C correspond to the benzopyran moiety, ring B to trihydroxyphenyl, and ring D to trihydroxybenzoate.

Methodology

Ligand parameters

The EGCG molecule was manually built using Gaussview (Millman & Keith, 2014). The structure was optimized using density functional theory with the m06–2× level (Zhao & Truhlar, 2007) and a 6–311++G (2d, 2p) basis set. We have previously used this combination of functional and basis set to optimize the structure of both organic and transition metal compounds with good approximation to experimental structures (Díaz-Gómez, Galindo-Murillo, & Cortés-Guzmán, 2017; Galindo-Murillo & Cheatham, 2014; Toledano-Magaña et al., 2017). Frequency calculations were performed and no negative values were found, ensuring the convergence of the structure. The resulting optimized structure was compared with experimental structures available from PDB entry codes 3OOB (X-ray, (Urusova et al., 2011)) and 2KDH (NMR, (Robertson, Li, & Sykes, 2009)). The RMSD using the benzopyran ring between the calculated structure and the X-ray structure was of .53 Å and .82 Å from the NMR data. Charge parameters for the optimized structure were generated using the R.E.D. server (Vanquelef et al., 2011). The general amber force field, GAFF (Wang, Wolf, Caldwell, Kollman, & Case, 2004), was used to describe the structural parameters.

Umbrella sampling

To study a possible intercalation mechanism of interaction, biased sampling was used by means of umbrella sampling to construct a pathway of intercalation between the DNA and rings A, B, and D of the EGCG ligand. This methodology consists of a series of intermediate steps along a pre-defined reaction coordinate that effectively samples an energy pathway with the aid of a biasing potential. Using a 36-mer with the sequence d(GCCTATAAACGCCTATAA), the EGCG molecule was manually intercalated between base-pairs A9-T28 and C10-G27 in three different configurations, mainly, the intercalations of ring A, B, and D as a starting structure. The separation restraint was centered on the heavy atoms of rings A, B, and D of the EGCG molecule and the center of mass of the two base-pairs forming the intercalation site. Multiple restrains were tested with different anchor points until we ensured that the ligand moved toward the major groove, while leaving the EGCG molecule to reorient freely once it was no longer in close proximity of the DNA chain. Because of the branching structure of the EGCG molecule, it is very unlikely that a similar process could take place through the minor groove, on which multiple steric clashes were observed when we tested a reaction coordinate through the minor groove. The sampling was carried out using .2 Å increment windows from 1 to 20 Å with a 10 kcal/mol Å harmonic restraint constant. Each window was sampled for 100 ns, and five independent copies of each window were performed to ensure statistical validity. The biasing potential was removed using the weighted-histogram analysis method, WHAM (Kumar, Rosenberg, Bouzida, Swendsen, & Kollman, 1992; Roux, 1995).

Simulation details

The 36-mer double strand DNA sequence used in this work has been previously studied (Lankas, Šponer, Langowski, & Cheatham, 2004). The DNA was built using the procedures previously described by the Ascona B-DNA Consortium, ABC (Beveridge, Cheatham, & Mezei, 2012; Beveridge et al., 2004; Dixit et al., 2004; Lavery et al., 2009; Pasi et al., 2014) and with the TIP3P (Jorgensen, Chandrasekhar, Madura, Impey, & Klein, 1983) water model using a truncated octahedral periodic cell and a 12 Å distance between the solute and the edge of the box. A total of 34 sodium ions were added to reach a neutral charge, and an excess of 48 sodium and chlorine ions was added to reach a final concentration of ~150 mM using the Joung-Cheatham set of ion parameters (Joung & Cheatham, 2008). After the model was built, the position of the ions was randomized by swapping a water molecule with an ion such that no ion was closer than 4 Å to another and the ions were farther than 6 Å away of the DNA–EGCG complex. This process was repeated four times to create four independent initial geometries, and each copy was run through the same minimization, equilibration, and production phase. The initial minimization process was performed using 1000 steps of steepest descent optimization switching to 1000 steps of conjugate gradients, applying a harmonic restraint of 25 kcal mol⁻¹ Å⁻¹ to the DNA and the EGCG molecules. Then, using a Langevin thermostat (Pastor, Brooks, & Szabo, 1988)with a collision frequency value of 1 ps⁻¹and the same restraints, the system was slowly heated from 100 to 300 K in 1000 ps of MD simulation. The restraint was decreased from 25 to .5 kcal mol⁻¹ Å⁻¹ in a series of equilibrium steps, each 2 ns long using constant pressure and temperature. The Berendsen barostat (pressure relaxation time set to 1 ps) was used (Berendsen, Postma, van Gunsteren, DiNola, & Haak, n.d.). A 2 fs time step was used for each equilibration step with a cutoff of 9.0 Å. MD simulations were performed using the CUDA version of pmemd available in AMBER 14(Case et al., 2014; Cheatham & Case, 2013; Götz et al., 2012; Le Grand, Götz, & Walker, 2013; Salomon-Ferrer et al., 2013). The PME method (Konerding, Cheatham, Kollman, & James, 1999) was used to calculate the long-range interaction and trajectory information was saved every 10 ps. Clustering analysis was performed using the hierarchical agglomerative algorithm as implemented in CPPTRAJ with a cutoff value of 3 Å between cluster centroids and using only heavy atoms.

EGCG–DNA binding energy analysis was calculated using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) methodology, (Hansson, Marelius, & Aqvist, 1998; Kollman et al., 2000) through the MMPBSA.py (Miller et al., 2012) application available in AmberTools16. A total of 5000 snapshots were used for the MM-PBSA analysis in each case.

Absorption spectra were calculated following a QM TD-DFT methodology with the hybrid functional of Adamo, PBE0(Adamo & Barone, 1999). This functional was chosen after running a series of tests calculation with other levels of theory (which included M062x, M11, B3PW91, cam-B3LYP, and PBE0) until we could find a good agreement with the reported experimental data. Additionally, this function has been tested thoroughly elsewhere (del Campo, Gázquez, Trickey, & Vela, 2012; Cárdenas-Jirón & Cortez, 2013). All QM calculations to obtain the UV-vis spectra were made using 6–311 + G (d, p). Implicit solvation was treated using the SMD dielectric method (Marenich, Cramer, & Truhlar, 2009). To reduce computation time, only residues 9–12, 25–28 and the EGCG ligand were used. This selected section of the DNA represents the full binding site of the EGCG molecule. Deleted bonds were capped with hydrogen atoms at both ends of the DNA chain. The remaining bases, sugars, and full backbone atoms were taken into account for the QM calculations. Quantum mechanical calculations were run using the G09-D.01 version of the Gaussian program (Frisch et al., 2009). Trajectory analysis was run using CPPTRAJ (Roe & Cheatham, 2013) as found in AmberTools 16.

Results and discussion

Intercalation pathway studies using umbrella sampling

The resulting free energy profile for each system and an average structure obtained from the minimum energy points is presented in Figure 2. The energy minimum for ring A is reached at 9.8 Å with a free energy value of −2.0 kcal/mol, 10.8 Å for ring B with a free energy value of −1.6 kcal/mol and 11 Å for ring D with a value of −1.2 kcal/mol. The intercalation process of ring A of EGCG decreases steadily as the benzopyran ring is removed from the intercalation site. At 4.8 Å there is a metastable intermediate with an average free energy value of 3.7 kcal/mol, which corresponds to a reorientation of the whole molecule that allows the opposing ring D to interact with the DNA backbone. The energetically stable (energy minima) configuration was found to be a partial stacking interaction for ring A, where the benzopyran moiety is still present between the ApC base-step. For rings B and D, the energy minima are with the corresponding ring outside the stacked bases. Free energy values for the stable interaction in each case show that ring A has preference over the other two conformations (even though with differences of .4 and .6 kcal/mol ,respectively, which is beyond the credibility of any molecular dynamics technique). The de-stacking interaction pathway for ring D also shows a metastable structure when the trihydroxybenzoate ring is at 5.7 Å with a value of −.5 kcal/mol. This structure corresponds to a stacking interaction where ring D is between the aromatic center of A9 and the geometric center between the base-pair of C10 and G27. When in this conformation, both ring A and B are forming hydrogen bonds with the oxygen atoms from the DNA backbone of the 5′-strand. The structural difference between the metastable and stable structures corresponds to ring D of the EGCG ligand shifting outside of the intercalation pocket and the nucleotides resetting to a canonical pairing conformation (Supplementary Figure S1). Visual inspection of the representative structures obtained from a clustering analysis for each case (Figure 2, bottom), suggested that ring A allowed a reconfiguration of rings B and D, facilitating the formation of non-bonded interactions with the DNA and hence, stabilizing the overall formation of the complex. As an initial test for this hypothesis, we measured the DNA bending angle, using the top, center, and bottom residues as anchor points (residue numbers 2, 35 for top anchor; 8, 9, 10, 28, 29, 30 for center and; 17, 20 for bottom anchor). The angle between the center of mass from these three anchors provides an overall measure of the amount of bending of the DNA. This analysis was applied to the trajectory (~10 000 frames) of the window with the energetically stable value, as mentioned before. The overall bending angle when the ring A is at the 9.8 Å position is 160.6° (±5.8); for ring B in the 10.8 Å position is 162.2° (±6.2), and for ring D in the 11 Å position is 155.3° (±10.6). The angle value for the metastable window for ring D located at 5.7 Å is 166.8 (±6.5).

Figure 2. Top: potential of mean force energy profiles for the intercalation process of rings A, B, and D, respectively. Bottom: most populated structure extracted using clustering analysis from 100 ns of sampling time from window 9.8 Å for ring A, 10.8 Å for ring B, and 11 Å for ring D.

Moreover, using the five independent trajectories corresponding to the lowest energy point in the PMF results, we calculated the number of contacts between the EGCG molecule and the DNA from the 100 ns trajectory information (which is the length of each independent window used for umbrella sampling). The average number of unique contacts using a threshold of 5.0 Å between the EGCG molecule and the DNA chain was 220 for ring A. Rings B and D present a broader population of number of contact between ~75–150 which agrees with our observation of ring A allowing the rest of the EGCG ligand to interact with the DNA.

The observation of a metastable partial intercalated state in the case of ring D at ~5 Å raised the question whether this configuration would remain stable in an extended unbiased molecular dynamics simulation. To explore this idea, we started with the same initial intercalated structures we used for the umbrella sampling analysis (each of the EGCG ligands intercalated between an ApCstep) and performed 10 unbiased independent simulations for 1.5 μs using the same equilibration and production protocols as described in the methods section (with randomized ions in each starting structure). To our surprise, the EGCG molecule intercalated with the trihydroxybenzoate moiety (ring D) remained stable in all of the 10 copies, whereas ring A and B detached from the intercalated site less than ~10 ns (Figure 3).

Figure 3. Distance values as measured from all the DNA residues to the center of mass of the EGCG ligand over a period of 1.5 µs. The starting structure was with the corresponding ring manually intercalated between base-pairs 9, 28, and 10, 27 (ApC step, the same structures used for the umbrella sampling simulations). Each line corresponds to one of 10 independent, unrestrained simulation.

When ring D is in the intercalated position it forms stacking interactions with both A9 and G27, the benzopyran (ring A) positions directly next to the DNA backbone forming hydrogen bonds with the phosphate oxygen atoms from residues 27 and 28, and ring B stabilizes with hydrogen bonds between the hydroxyl oxygen atoms and N6 of adenine 6 (Figure 4). The sum of the previously described interactions is possible when ring D is in the intercalated position mainly because the branched structure of the EGCG ligand places both rings A and B in an optimum configuration to interact with the DNA backbone and the base-pairs in the interaction site. When ring A is stacked with the flanking bases, the aromatic ring of the benzopyran forms a stacking interaction with adenine 9 (A9) and guanine 27 (G27) in a similar manner as with the ring D position, although, since the benzopyran is a dual ring structure, rings B and D are facing away from the DNA, which do not allow the formation of stabilizing non-bonded interactions. Because rings B and D are freely rotating aromatic rings, the hindrance of these fluctuations causes ring A to slide away from the intercalation site. In the case of ring B, and in a similar way as ring D, the aromatic ring forms stacking interactions with A9 and G27. The problem in this case is that the aromatic ring B is most stabilized the closer it can get in contacts with the aromatic rings of both adenine and guanine present in the intercalation pocket. This means that the ligand is trying to push itself deeper in the intercalation pocket, but the branching structure of the ligand forms steric clashes that do not allow a proper π–π overlap interaction to occur between adenine – trihydroxyphenyl – guanine. Eventually, the fluctuations of the system overcome the weak formation of the complex and the ligand detaches from the intercalation site. This does not happen with ring D in the intercalated mode since the ketone group present in the benzoate structure allows the aromatic ring to be inserted deeper in the pocket and form stable stacking interactions before steric clashes occur by rings A and B.

Figure 4. View from 5′ → 3′ (left) and 3′ → 5′ (right) of the representative structure from the top populated cluster extracted using the 10 independent simulations of the EGCG ligand with the ring D in the intercalation position.

The binding free energy for the three complexes from the umbrella sampling simulations and the 10 independent unbiased MD simulations of ring D (Figure 3) is presented in Table 1. For each of the umbrella sampling simulations, the binding energy was calculated using a representative structure obtained from the lowest energy window by means of a clustering algorithm as described in the methodology section. Similarly, the binding energy of the unbiased simulations was calculated from a representative structure extracted from the clustering analysis using the 10 independent set of trajectories. The binding energy obtained from the unbiased simulation present the lowest value with a difference of ~20 kcal/mol from the structures obtained from the umbrella sampling studies.

Table 1. Binding energy values between DNA in the four proposed binding configurations and the maximum absorption wavelength as reported in the literature for EGCG and calculated using TD-DFT. The binding energy for the unbiased ‘Complex from MD – ring D’ is the average of the 10 independent simulations. Numbers in parenthesis represent standard deviation.

System	MM-PBSA (kcal/mol)	λmax (nm)
EGCG experimental value		272–275
EGCG calculated value		283
Complex from MD – ring D (unbiased)	−23.57 (±1.7)	308
Umbrella sampling – ring A structure	−4.9 (±5.1)	274
Umbrella sampling – ring B structure	−5.4 (±3.0)	274
Umbrella sampling – ring D structure	−1.8 (±2.6)	293

The energy values obtained from the umbrella sampling simulations are within ~4 kcal/mol within each other, which is expected since their binding modes are similar: In each case, the EGCG molecule is interacting within the major groove of the DNA with both the sugar moiety and the phosphate groups, and little interaction with the actual DNA bases is present. No real binding difference could be proposed based on the umbrella sampling results. Contrary to this, the EGCG–DNA binding mode of the unbiased simulations with the ring D in intercalation presents a binding energy value which is more ~15 kcal/mol stronger than the binding poses observed from the umbrella sampling simulations and is within the range of the computed free binding energy of similar sized molecules interacting within the grooves of DNA (Díaz-Gómez et al., 2017; Galindo-Murillo & Cheatham, 2014; Galindo-Murillo et al., 2015; García-Ramos, Galindo-Murillo, Cortés-Guzmán, & Ruiz-Azuara, 2013). It is important to mention that the presented energy values are over-estimates not only due to the PBSA which tends to over-estimate binding free energies, but also the loss of translational and rotational entropy (~12 kcal/mol) were not taken into account in the free energy calculations (Homeyer & Gohlke, 2012; Štefl et al., 2003).

Spectroscopic calculations

To date, no published X-ray or NMR structure of the interaction of any type of catechin molecule with nucleic acids is available, although there is strong experimental evidence that these molecules do bind with DNA and RNA (Kuzuhara et al., 2006). Dasgupta and collaborators reported the experimental absorption profile for EGCG with increasing concentrations of calf-thymus DNA (Supplementary Figure 4 of reference (Ghosh et al., 2008)). This experimental titration of the EGCG molecule shows the hypochromic effect as DNA concentration increases, and bathochromic shift of the λ_max by ~10 nm. The absorption spectra were computed at the QM TD-DFT level for each of the three different conformation studied using the umbrella sampling methodology. The structures used in the QM calculations are the ones shown in Figure 2 (bottom), which are extracted after a clustering analysis from the lowest energy window. Additionally, we also calculated the spectra for the most populated structure obtained in the unbiased MD simulation (with ring D in intercalation mode). Comparison between the experimental and calculated UV-vis absorption information is presented in Figure 4. The calculated absorption’s maxima obtained for the EGCG molecule is at 279 nm, which is off by 5 nm from the reported experimental values at ~272–275 nm (in water). Ring A and ring B structures both show a λ_max = 274 (Table 1) with a reduction in the absorbance of ~.005. In contrast, ring D (from US) and the intercalated complex showed hypochromic effect of ~.04 and redshift of 19 and 34 nm, respectively. The shifts of λ_max and absorption from the complex structure obtained from the MD trajectories are observed in the experimental EGCG–DNA titration experiments, suggesting that the experimental binding mode is with the trihydroxybenzoate intercalating between the DNA base-pairs (Figure 5).

Figure 5. Normalized absorbance spectra for the EGCG molecule (circle). DNA–EGCG complex structure (diamonds) using the representative structure from a clustering analysis from the simulation where ring D is in the intercalated position within the DNA (Figure 4) and the structure from the PMF, ring D (cross ticks). Spectra for the PMF’s from ring A (up triangle tick) and ring B (down triangle tick) are extracted from a clustering analysis using the lowest energy window. Absorption spectra were calculated using the PBE0 functional with the SMD implicit solvation model in water. Experimental values of the EGCG molecule are shown as lines from references: (19, 20).

Sequence dependence analysis

The stable intercalation configuration from the trihydroxybenzoate (ring D) moiety of the EGCG molecule at the ApC base-step opened a new set of ideas regarding possible DNA sequence selectivity. To test this, we used the most populated structure of our previous simulation presented in Figure 4 as starting structure and manually switched the top and bottom base-pairs of the intercalation site with every possible combination of base-steps. The performed simulations are shown in Table 2.The resulting structures were minimized following the same procedure as described in the methods section. Three independent simulations for each base-step were sampled for 5 μs. The DNA–EGCG major groove intercalation complex was stable for the entire sampling time-tested, regardless of the sequence and overall, structural deformation of the DNA was only detected to a significant amount as fraying events on both ends of the DNA which are routinely observed for double-stranded DNA simulations (Galindo-Murillo, Roe, & Cheatham, 2014a; Galindo-Murillo, Roe, & Cheatham, 2014b; Zgarbová, Otyepka, Šponer, Lankaš, & Jurečka, 2014). A normalized histogram of the RMSD for the three aggregated trajectories of each system is shown in Figure 6, top. The simulations with the base-steps ApA, ApG, ApT, CpA, and GpG show a stable intercalation complex populated between ~2–3 Å. The additional structures that are being populated for CpC, GpA, GpT between ~2.5–4 Å and to a lesser extent with TpT and CpG are due to long-lived frayed event deformations present at both ends of the terminal base-pairs. This can be confirmed by means of a root mean square fluctuations (RMSF) analysis on each residue (Figure 6 bottom). In the figure, we notice that residues 5, 6, 15, and 16 show fluctuations between ~.1–.5 Å for all the tested systems, whereas residues present at both ends of the DNA chain show fluctuations in the ~4–8 Å range. This observation also indicates a stable complex is present regardless the sequence, and the different populations present in the RMSD analysis are not involved in the DNA–EGCG interaction.

Table 2. Simulations using different base-steps in the intercalation site. The original sequence corresponds to d(GCCTATAA-X-I-Y-GCCTATAA), where the X and Y positions are the swapped base-pairs and I marks the intercalation site. The sequence column describes the switched base-step.

Sequence	Binding energy (kcal/mol)
ApA	−22.15 ± .12
ApG	−23.40 ± .36
ApT	.23.02 ± .64
CpA	−21.88 ± .61
CpC	−22.21 ± .47
CpG	−23.84 ± .87
GpA	−21.84 ± .16
GpG	−24.11 ± .33
GpT	−22.97 ± .39
TpT	−21.26 ± .06

Figure 6. Top: normalized histogram of the RMSD values using the first frame as reference over the three aggregated independent copies for each system. Black line includes all residues; red line is RMSD without the terminal base-pairs at each end of the DNA. Bottom: RMSF by residue of each system using the full aggregated trajectory.

Energy analysis between the EGCG and the DNA using the MM-PBSA protocol and different sequences does not provide a clear trend (Table 2). The results all indicate an energetically favorable complex as observed in the RMSD and RMSF data. Among the data-set, the binding energy values of stable structure as suggested by the binding values are observed with the GpG, CpG, and ApG present the lowest energy value and TpT, GpA, and CpA the highest values, although, the average difference between these six simulations is ~2.8 kcal/mol. Hydrogen bonding analysis using the 10 simulations showed that interactions with the backbone provide the main stabilization point for the complex. The interaction between the hydroxyl group of ring A (position 5, Figure 1) with the phosphate group of residue 15 is present in all simulations (Supplementary Table S2). Additionally, the hydroxyl groups of ring B also form hydrogen bonding with the phosphate group on the opposite side located in residue 5. The most populated structure for each individual system using the full aggregated trajectory shows the EGCG ligand in the intercalation pocket in each case (Supplementary information S2).

Visual inspection of the EGCG molecule in the intercalation site indicates that the ligand does not show significant conformational changes throughout the simulations regardless of the sequence. This observation was further tested by measuring the flexible dihedral angles between the trihydroxyphenyl, trihydroxybenzoate, and the benzopyran rings (the bonds between ring B to C and D to C, Supplementary Figure S3). The dihedral angle between the trihydroxyphenyl moiety and the benzopyran structure is populated in a balanced range between 101° and 270° for the entire aggregated simulations that corresponds to the expected free rotation of the aromatic ring. For the trihydroxybenzoate ring, two dihedrals were monitored located in the methyl ester bond connecting the aromatic trihydroxyphenyl ring with the benzopyran (Supplementary Figure S3). Both dihedral bonds populate values close to ~167° and ~187°, respectively, which correspond to the conformation previously observed (Figure 4) and no major shifts are observed confirming a lack of conformational shifts of the EGCG ligand.

Conclusions

The branched structure of the EGCG molecule limits the way it can interact with double-stranded DNA. Numerous steric clashes make interaction of the ligand through the minor groove of the DNA energetically unfavorable hence we hypothesized that the interaction would take place through the major groove. We have then proposed three possible binding configurations based on the number of aromatic rings present in the ligand. Using enhanced sampling molecular dynamics, we explored the possibility of an intercalation interaction between each of the aromatic sites of the EGCG ligand between an ApC DNA step. The results showed energetically stable configurations outside the intercalation pockets in all the cases, however, the EGCG–DNA interaction using the trihydroxybenzoate (ring D) showed an additional stable configuration at a distance of ~5 Å between the center of mass of both EGCG ligand and DNA axis. Upon visual examination of this structure, we confirmed that the ligand was in full intercalation mode between the DNA base-pairs. Unbiased, unrestrained simulations of the DNA-ligand complex starting from the manual intercalation of each of the EGCG aromatic rings yielded a stable complex structure of the trihydroxybenzoate moiety with DNA after 10 independent trajectories (each > 1.5 μs in length), whereas intercalation structures using the trihydroxyphenyl or the benzopyran moiety were unstable and the intercalation interaction was lost in ~50–100 ns of simulation time in all of the 10 independent runs.

It is important to mention that the stable intercalation binding mode observed between the trihydroxybenzoate and the DNA could potentially be biased by the known over-stacking effect present in the parm99-bsc0 AMBER force field as reported by Garcia and collaborators (Chen & Garcia, 2013). A follow-up and more recent work performed by Elcock and collaborators studying the stacking energies of dinucleotide monophosphates mention the possibility of a more favorable for RNA than DNA, although they conclude that this observation is dependent of the presence of T/U bases and the experimental databases on which comparison is made. We hypothesize that even though the over estimation of stacking interactions is in fact, present, and potentially facilitating the intercalation of the EGCG with DNA, this biasing was not observed on two of the three binding poses we explored. In fact, when manually intercalating rings A and B of the EGCG molecule, the intercalation pose was lost in a ~200 ns time scale, even though, 10 independent simulations were run for each of the systems. Further research in this particular stacking subject will be addressed in future work, noting that in more recent work on dinucleotides by our lab with the various AMBER force fields (work under review), the over-stacking observed by Elcock is less profound than previously reported.

The QM calculation of absorption spectra of the stable intercalation complex showed a hypochromic effect with a redshift of ~10 nm, which is consistent with what is observed in experiments of the EGCG molecule in the presence of calf-thymus absorption experiments. Possible DNA sequence preference by the ligand was studied, although no clear trend was observed. Our observations have generated atomistic details about the major groove intercalation process of a known antioxidant molecule, moving forward the studies of how DNA interacts with molecules.

Disclosure statement

No potential conflict of interest was reported by the authors.

Funding

This work was supported by the Blue Waters sustained-petascale computing project of National Science Foundation [grant number NSF OCI 07-25070], [grant number PRAC OCI-1036208]; the NSF Extreme Science and Engineering Discovery Environment [grant number XSEDE, OCI-1053575], allocation [grant number MCA01S027P]; the Center for High Performance Computing at the University of Utah.

Supplementary material

The supplementary material for this paper is available online at https://doi.org/10.1080/07391102.2017.1389306.

Acknowledgments

We would like to thank James Robertson for his comments about the manuscript.