The binding modes of cationic porphyrin-anthraquinone hybrids to DNA duplexes: in silico study

Cationic porphyrin-anthraquinone hybrids bearing peripheral substituents, either pyridine, imidazole, or pyrazole rings have been investigated for their binding mode to DNA duplexes. The four kinds of DNA duplexes were used, which represent intercalation and groove binding modes. AutoDock 4.2 was used to dock nine hybrid compounds to four DNA duplexes, while monitoring of conformational changes of four best hybrid–DNA complexes during 2 ns was performed by Amber9 molecular dynamics package. The binding energy calculation of best four complexes was then carried out using MMPBSA method. The hybrid compounds interacted to DNA duplexes through intercalation and groove binding modes. The minor groove binding of DNA was energetically preferred by cationic porphyrin hybrids due to favorable electrostatic and van der Waals interactions. Both electrostatic and van der Waals contributions were able to distinguish the binding mode of porphyrin hybrid to DNA duplexes.


Introduction
Cationic porphyrins are known to have strong interaction with DNA, high photonuclease and photodynamic therapy activities, and high inhibitory activity of telomerase (Fiel, 1989;Guliaev & Leontis, 1999;Kim, Gleason-Guzman, Izbicka, Nishioka, & Hurley, 2003;Vedachalam et al., 2011;Wang, Zhang, Wu, & Li, 2007;Yamamoto, Tjahjono, Yoshioka, & Inoue, 2003). Those favorable properties are greatly encouraged by the tight binding of porphyrin to DNA. Essentially, cationic porphyrins bind to DNA by three distinct modes: intercalation, external binding, and external binding with self-stacking. A number of factors are known to influence the mode and sequence selectivity of DNA binding such as ionic strength, the presence and nature of charges of metal complexed in porphyrin core, and the nature and position of substituents on the porphyrin ring system (Sari, Battioni, Dupre, Mansuy, & Le Pecq, 1990;Wu et al., 2006).
Introduction of several functional groups such as pyrimidine, imidazole, and pyrazole in the periphery of porphyrin core has been shown to increase porphyrin binding to DNA. In our previous work, cationic porphyrin bearing imidazole or pyridine substituents binds outside to the minor groove of calf thymus DNA (ctDNA) while those with pyrazole substituent intercalate into ctDNA. The fitting of pyridinium ring into the DNA base pairs causes larger entropy changes than those of the pyrazolium ring. It was implied that structural differences in meso-substituents greatly vary thermodynamic parameters of the binding of porphyrin-DNA complex (Tjahjono, Akutsu, Yoshioka, & Inoue, 1999).
In the meantime, in the effort of increasing porphyrin binding to DNA, a search for particular chemical units to be introduced in the meso-20-position have been continuously increased. Among them, quinoxaline, aminoquinoline, alkaloid, peptide, acridine, phenylpiperazine, and anthraquinone groups are reportedly successful to be conjugated on the porphyrins and showing interesting binding properties to DNA (Ishikawa, Yamashita, & Uno, 2001;Jia, Jiang, Wang, & Li, 2006;Kumar et al., 2013;Li, Fedorova, Trumble, Fletcher, Czuchajowski, 1997;Maraval et al., 2003;Orosz et al., 2013;Veverková, Záruba, Koukolová, & Král, 2010;Zhao et al., 2008). In the study of cationic porphyrinanthraquinone hybrids, Zhao et al. (2008) implied that there was intense interaction between the anthraquinone group conjugated to porphyrin core and DNA. While cationic porphyrin moiety binds to DNA in outside mode due to its excellent planarity and the hindrance, they proposed that anthraquinone plane binds to DNA in partial intercalative mode. They also suggest that both porphyrin and anthraquinone in the hybrid interact with DNA rather strongly (Zhao et al., 2008).
Despite of those facts, little is known regarding the mechanisms of the binding of cationic porphyrin-anthraquinone hybrids to DNA in the molecular level. The present study is conducted to understand DNA binding of the porphyrin hybrid bearing either pyridine, imidazole, or pyrazole groups at the meso-positions by using docking and molecular dynamics approaches. While the use of molecular docking has been evidenced to be capable of reproducing crystal structure of minor groove and intercalation agents of DNA complexes, molecular dynamics is widely used in monitoring conformational changes in ligand-macromolecule complex during a certain period of time (Holt, Chaires, & Trent, 2008). Using molecular dynamics trajectory, prediction of binding energy of a complex is widely accepted by MMPBSA method. Thus, the present study explains the mechanism of action and thermodynamics properties influencing dynamics of cationic porphyrins-anthraquinone hybrid towards DNA duplexes either through intercalation, minor groove, or a combination of both.

Macromolecule preparation
Crystal structures of the starting DNA receptor were downloaded from the Protein Data Bank with pdb code 2DND and 1D64 for dodecamers d(CGCAAATTTGCG) 2 and d(CGCGAATTCGCG) 2 , herein denoted 2DND and 1D64, respectively, which represent the mode of minor groove binding of DNA duplex, and the hexamers d(CGATCG) 2 pdb code 1Z3F and d(TGATCA) 2 pdb code 182D, herein denoted 1Z3F and 182D, respectively, which represent the mode of intercalation binding of duplex DNA. At each macromolecular target, the natural ligand and water molecules were removed, and polar hydrogen and Kollman charges were assigned. Figure 1 shows nine porphyrin hybrids used in the present study, i.e. mono-H 2 ImP-AQ, mono-H 2 PyP-AQ, mono-H 2 PzP-AQ, bis-H 2 ImP-AQ, bis-H 2 PyP-AQ, bis-H 2 PzP-AQ, tris-H 2 ImP-AQ, tris-H 2 PyP-AQ, and tris-H 2 PzP-AQ. Each ligand structure was built with the GaussView. Geometry optimization was performed using Gaussian03 on 6-31G* basis set using the Becke threeparameter Lee-Yang-Parr (B3LYP) hybrid density functional theory. This level of theory is suitable for such tasks (Brovarets, Yurenko, Dubey, & Hovorun, 2012;Brovarets, Yurenko, & Hovorun, 2013;Ponomareva, Yurenko, Zhurakivsky, van Mourik, & Hovorun, 2013). Structure with the lowest energy conformation was selected for further docking study. Each ligand molecule was then prepared with the program, AutoDock Tools with maximum torsion (AutoDockTools, 2007).

Molecular docking
AutoDock 4.2 was used to dock nine ligands in four different DNA duplexes, while AutoDock Tools 1.5.6 was employed to establish the AutoGrid points of porphyrin hybrid-DNA complexes (Goodsell, Morris, & Olson, 1996;Morris et al., 1998). The grid maps were established by centering the grid box on either the minor groove or the intercalation site and comprised 62 × 62 × 62 points of 0.375 Å spacing. The Lamarckian genetic algorithm was used with population size of individuals = 150, maximum number of energy evaluations = 2.5 × 10 6 , maximum number of generations = 27,000, elitism value = 1, mutation rate = 0.02, and crossover rate = 0.80. Docking runs of 50 were employed in all the calculations with step sizes of 0.2 Å for translations and with orientations and torsions of 5.0°. All other docking parameters were left at the default values.

Molecular dynamics simulation
Molecular dynamics simulations were performed with the SANDER module of Amber v9.0 program (Case et al., 2006) on four porphyrin hybrid-DNA complexes, which were best docking conformation of the ligand in each DNA receptor. AMBER force field ff99 and Gaff (Yang et al., 2006) were used to parameterize the DNA and ligand, respectively, with the assign of partial atomic charges AM1-BCC for ligand. The addition of Na + ions in the "xleap" done to neutralize system. Solvation was done with truncated octahedron TIP3P box and the distance of 10 Å from the edge of complex atoms. Energy minimization and MD simulations were carried out as follows. First, energy minimization was carried out using 1500 steps of steepest descents and 1500 steps of conjugate gradients with the DNA and ligand restrained (force constant 500 kcal/mol Å 2 ). Next, the system was left free to move and the complex was minimized by 1500 steps of steepest descents and 1500 steps of conjugate gradients. The system was heated from 0 to 300 K for 50 ps with a time step of 0.0005 ps and positional restraints for the DNA and ligand atoms with a 10 kcal/mol Å 2 force constant and constant volume molecular dynamics. Berendsen (Berendsen, Postma, van Gunsteren, DiNola, & Haak, 1984) coupling of 0.2 ps was used to maintain the temperature and system pressure. Bond lengths involving hydrogen atoms constrained with SHAKE (Ryckaert, Ciccotti, & Berendsen, 1977) algorithm with 2 fs integration time step. PME method (Darden, York, & Pedersen, 1993) was used for the treatment of electrostatic interactions with a cut-off distance 8.0 Å for nonbonded interactions. Then, the temperature was maintained at 300 K during equilibration and production for 2 ns using periodic boundary condition and constant number of particles, pressure, and temperature ensemble. The reference pressure and temperature were set to 1 atm and 300 K, respectively. Verification of system stability analysis was performed by the convergence of energy, temperature, pressure, and root mean square deviation of the DNA and the ligand (RMSD). Visualization was done with VMD program 1.9 (Humphrey, Dalke, & Schulten, 1996).

MM/PBSA calculations
According to MMPBSA method (Kollman et al., 2000), the binding free energy of each system is formulated with the following equation: where G complex , G rec , and G ligand is a complex, receptor, and ligand free energy, respectively. ΔE MM is the interaction energy between the receptor and the ligand donated by the bonded terms (bond energy, angle energy, and torsion energy), and nonbonded terms (van der Waals energy and electrostatic energy). ΔG PB and ΔG SA are the polar and nonpolar contributions to desolvation upon ligand binding, respectively; and −TΔS is the conformation entropic contribution at temperature T. Polar desolvation energies were calculated by solving the Poisson-Boltzmann (PB) equation with a grid size of 0.5 Å. The dielectric permittivity of solute interior was set to 1 which was consistent to the nonpolarizable force field used, while the dielectric permittivity of 80 was used for describing water (Hou et al., 2010;Špačková et al., 2003). The nonpolar contribution was calculated by solvent accessible surface area with solvent probe radius of 1.4 Å. Free energy calculations were performed with the python modules as implemented in Amber9. The normal mode analysis (Kollman et al., 2000) was carried out to calculate the entropy change of conformational ligand binding (−TΔS). Free energy of DNAligand binding was calculated based on 200 snapshots taken from 1 to 2 ns MD simulation trajectories of the complex.
Results and discussion Figure 2 shows energy-minimized geometries of tris-H 2 ImP-AQ and tris-H 2 PyP-AQ. In each optimized structure, the anthraquinone group was perpendicular to the porphyrin macrocycle. The long linear chain connecting porphyrin and anthraquinone ends allows the two moieties to contribute simultaneously upon interaction with DNA duplexes. The similar geometries were also observed for other ligands ( Figure S1).

Binding modes of porphyrin hybrid-DNA duplex complexes
Each ligand was docked on each of the four duplex DNAs. Docking conformations with the lowest energy of the best cluster are taken for next analysis. The docking results of nine ligands on the four DNA duplexes are shown in Table 1. The nine ligands can bind both to the minor groove of dodecamers and can be intercalated to the hexamer. Table 1 shows that bis-H 2 PyP-AQ is best docked on minor groove binding representing DNAs of both 2DND and 1D64 with binding energies of −14.84 and −11.54 kcal/mol, respectively. In both conformations, the crescent shape of the porphyrin hybrid complements the bend of the DNA double helix with close hydrophobic interaction (Martinez and Iverson, 2012). As for 2DND, bis-H 2 PyP-AQ also bind through minor groove of 1D64 with lowest binding energy; however, when bound to 1D64, the conformation of the formed complex was not as good as 2DND complex with binding energy difference as many as 3.3 kcal/mol. It was evidenced that porphyrin hybrid binding through minor groove of duplex DNA prefers AT-rich sequence (Pasternack, Gibbs, & Villafranca, 1983). Interestingly, the ligand was also docked well on intercalation binding representing DNA i.e. 1Z3F with binding energy −10.52 kcal/mol, only slightly different from the binding energy of mono-H 2 PzP-AQ which has the lowest binding energy of −10.66 kcal/mol on the hexamer. The ring system of porphyrin core stacks well with base pairs in the intercalation complex involving strong π − π interaction between relatively hydrophobic and planar porphyrin macrocycle with adjacent groups of DNA chain while anthraquinone end tends to interact through minor groove of duplex DNA (Figure 3). Thus, bis-H 2 PyP-AQ has potential to be both as minor groove and intercalation binder of the DNA duplex. The π − π interaction was represented by yellow color in Figure S2.
The binding energies of porphyrin hybrid-DNA complexes formed by binding through minor groove were much more negative compared to that formed by binding through intercalation mode. It might be due to a lack of planarity and steric hindrance in the periphery of the planar porphyrin region that inhibits intercalation, thereby the binding of porphyrin hybrids favors minor groove of DNA. This suggests that the mode of porphyrin hybrids binding through the minor groove of DNA duplex is more energetically favored than the ligand binding via intercalation mode. Additionally, ligands bearing one (mono-) and two (bis-) substituents have relatively lower binding energies than the ligand binding energy of those with three substituents (tris-), no matter on which duplex DNA the ligand was bound. Again it appears that sterical hindrance due to the side groups on the tris compounds induced constraint to the ligands to attain the best geometry upon interaction with duplex DNA. This was confirmed by the pattern of mono-and bis-substituents ligand interactions in which the nonbonded interactions were contributed by groups other than the side groups of porphyrin macrocycle. In general, the nine ligands prefer interaction with 2DND than with 1D64. It seems that there has been a more intense interaction for minor groove binding between the anthraquinone group of porphyrin hybrids with thymine residues in dodecamer 2DND compared to interaction of anthraquinone group of hybrid with cytosin residue of 1D64.
Moreover, when the hybrids were docked to intercalation binding representing DNAs, the binding energy of porphyrin hybrids to 182D was more negative (stronger affinity) than those to 1Z3F. This suggests that intercalation of porphyrin hybrids in these two DNA duplexes prefers AT-rich sequence of DNA. In general, ligand conformation when bound to 182D left anthraquinone group as intercalator, however when bound to 1Z3F, the ligand left porphyrin macrocycle as intercalator. Figure 4 shows the conformation of mono-H 2 PyP-AQ intercalated into base pairs of d(CGATCG) 2 and base pairs of d (TGATCA) 2 . The π − π interaction was represented by yellow color in Figure S3.
Macrocycle porphyrin of tris-H 2 PyP-AQ in the interaction with dodecamer 2DND leads out while  porphyrin ring of tris-H 2 ImP-AQ as well as anthraquinone group of tris-H 2 PzP-AQ both leading to the major groove of 1D64 dan hexamer 1Z3F, respectively. The existence of both porphyrin and anthraquinone group that leads to the major groove of DNA duplex seems to be a energetically unfavorable, characterized by more positive binding energy as observed in the tris-H 2 PzP-AQ, and tris-H 2 PyP-AQ to 1D64 ( Figure 5). Interactions in the major groove of hexamer d(CGATCG) 2 (pdb code 1Z3F) were also observed in almost all the ligands. Furthermore, we also observed patterns of molecular dynamics on the most stable of ligand-DNA complexes, i.e. bis-H 2 PyP-AQ and 2DND, bis-H 2 PyP-AQ and 1D64, mono-H 2 PzP-AQ and 1Z3F, as well as bis-H 2 PyP-AQ and 182D. Figures 6 and 7 show the total energy profiles of each complex with respect to time during 2 ns.
In each complex, the average total energy is relatively stable during the 2 ns, but there are striking differences in the four complexes. In the first two complexes, which are the complexes formed through minor groove binding, a total energy of the two complexes were nearly equal, while the total energy in the third complex formed via the intercalation binding mode is much more positive than that in both complexes formed through the minor groove binding. This suggests that the mode of ligand binding through the minor groove of  DNA duplex (two first complexes) is more energetically favored than the ligand binding via intercalation mode (third complex). In addition, we also observed RMSD values of each complex formed as shown in Figure 8. The four complexes underwent a minor movement, which has quite reasonable RMSD values with an average below 3 Å. Table 2 shows the predicted binding affinities and the energy components of the four porphyrin hybrid-DNA complexes. Clearly, the predicted binding free energy of the bis-H 2 PyP-AQ to 1D64 complex (−78.17 kcal/mol) was more favorable than that of the bis-H 2 PyP-AQ to 2DND complex (−69.48 kcal/mol). Additionally, the binding free energy of the mono-H 2 PzP-AQ to 1Z3F (−52.14 kcal/mol) was less favorable than that of the aforementioned two first complexes, showing that intercalation of mono-H 2 PzP-AQ to 1Z3F was less favorable compared to minor groove binding of the bis-H 2 PyP-AQ to both 2DND and 1D64.

Binding free energies predicted by MM/PBSA
Furthermore, we examined the four individual energy components (ΔE ele , ΔE vdw , ΔG PB , and ΔG SA ) in order to    acquire the role of each energy term on the binding energies. In Table 2, the binding bis-H 2 PyP-AQ to 2DND was largely determined by both the electrostatic (ΔE ele ) and the van der Waals (ΔE vdw ) energy terms. The electrostatic forces which involved in porphyrin-DNA interactions are mainly contributed by positively charged meso-substituents of porphyrin and negatively charged phosphate oxygens of DNA (Tjahjono et al., 2006). The electrostatic interaction (ΔE ele ) of bis-H 2 PyP-AQ to 2DND complex (−1048.69 kcal/mol) was slightly lower than that of bis-H 2 PyP-AQ to 1D64; however, it was considerably lower compared to those of intercalationrepresenting complexes (−279.19 and −597.69 kcal/mol for mono-H 2 PzP-AQ to 1Z3F and bis-H 2 ImP-AQ to 182D, respectively). It was also the case in the van der Waals contribution (ΔE vdw ) of the bis-H 2 PyP-AQ and 2DND complex (−83.97 kcal/mol), which was lower than that of bis-H 2 PyP-AQ and 1D64 complex (−78.77 kcal/mol). However, the van der Waals term increases sharply in intercalation-representing complexes (−69.88 and −68.81 kcal/mol for mono-H 2 PzP-AQ to 1Z3F and bis-H 2 ImP-AQ to 182D, respectively). As such, it is clear that the van der Waals contribution between the two binding modes was distinctly different. From this point of view, ligand binding through minor groove is favorable than that through intercalation, and both electrostatic and van der waals term are essential to distinguish between minor groove and intercalation binding modes.
The electrostatic term obtained by both electrostatic and polar contribution of desolvation (ΔE ele + ΔG PB ) of the porphyrin hybrid-DNA complexes were unfavorable for the binding of bis-H 2 PyP-AQ (30.48 kcal/mol and 17.96 kcal/mol to 2DND and 1D64, respectively, and 29.35 and 23.27 kcal/mol for mono-H 2 PzP-AQ to 1Z3F and bis-H 2 ImP-AQ to 182D, respectively. The electrostatic term contributes more positive to the total enthalpy of bis-H 2 PyP-AQ to 2DND by 12.52 kcal/mol compared to that of bis-H 2 PyP-AQ to 1D64. The nonpolar desolvation energy term (ΔG SA ) was different for the two modes of interaction among the porphyrin hybrid-DNA complexes (−44.90 and −42.92 kcal/mol for minor groove recognition of bis-H 2 PyP-AQ to 2DND and bis-H 2 PyP-AQ to 1D64, respectively, and −36.71 and −37.62 kcal/mol for intercalation of mono-H 2 PzP-AQ to 1Z3F and bis-H 2 ImP-AQ to 182D, respectively). Furthermore, the van der Waals contribution difference between bis-H 2 PyP-AQ to 2DND and bis-H 2 PyP-AQ to 1D64 complexes was only 5.20 kcal/mol; however, the difference of that among the bis-H 2 PyP-AQ to 2DND with mono-H 2 PzP-AQ to 1Z3F and bis-H 2 ImP-AQ to 182D complexes were 14.09, and 15.16 kcal/mol, respectively. These results indicate that the preference of the binding mode between minor groove and intercalation was significantly determined by the van der Waals contribution.

Conclusions
A series of cationic porphyrin-anthraquinone hybrids are able to bind to DNA duplex through minor groove, intercalation, as well as major groove binding. Compared to other binding modes, minor groove binding is apparently favorable energetically. In term of ligand structure, DNA binding prefers side ring with less number of substituents in order to dismiss sterical hindrance upon ligand-DNA interaction. Molecular dynamics studies suggest that the porphyrin hybrid-DNA duplex complexes undergo only a minor deviation during 2 ns. It was found that both electrostatic and van der Waals contribution were crucial in the binding mode preference of porphyrin hybrid to DNA duplexes.