Molecular Modeling Studies of Interaction Between Plasmid DNA (pBR322) and Dendritic Antioxidants

It is known that excessively produced reactive oxygen species (ROS) like free radicals as a consequence of normal metabolism and by inflammatory mediators during (chronic) inflammation can attack cellular materials like DNA, lipid, protein, membranes and genes, and cause oxidative stress to the cells. In the case of DNA, oxidative damages caused by reactive free radicals include single or double strand breaks, base modifications, or rearrangements, which result in alteration of gene expression, cell growth and differentiation leading to mutagenesis as well as carcinogenesis [2]. Therefore, reduction of such damages by efficient management of free radicals can result in a reduced risk of diseases including cancer.


Introduction
It is known that excessively produced reactive oxygen species (ROS) like free radicals as a consequence of normal metabolism and by inflammatory mediators during (chronic) inflammation can attack cellular materials like DNA, lipid, protein, membranes and genes, and cause oxidative stress to the cells. In the case of DNA, oxidative damages caused by reactive free radicals include single or double strand breaks, base modifications, or rearrangements, which result in alteration of gene expression, cell growth and differentiation leading to mutagenesis as well as carcinogenesis [2]. Therefore, reduction of such damages by efficient management of free radicals can result in a reduced risk of diseases including cancer.
Antioxidants are compounds, which donate electron(s) or hydrogen atom(s) to reactive free radicals, thus inhibiting chain reactions and preventing oxidation of other molecules. Although the mechanism through which antioxidant works have not been thoroughly understood yet but significant advances have been made in understanding how they contribute to reduction of oxidative stress and prevent the pathogenesis of various diseases and aging. Our system consists of compartments with different hydrophilicity/ hydrophobicity. Therefore, the distribution of antioxidants in those compartments normally depends on their water-lipid solubility, which is dictated by their chemical structure. The chemical structures affect not only the local concentrations of antioxidant but also interactions with particular biomolecules at the site [3,4]. Thus, the capacity of an antioxidant to interact well with a biomolecule might be important to produce potent antioxidant activity towards the biomolecule and to regulate possibly other cellular events. There are reports stating that potent plant-based phenolic antioxidants like flavonoids bind well with DNA. Among the flavonoids, quercetin is one of the most widely studied flavonoids and currently being exploited as a chemo-preventive agent. Quercetin is known to bind to DNA through intercalation and external binding, and prevent DNA from free radical damage [5]. For quercetin, the intercalation was conceivably facilitated by its planar molecular conformation [6][7][8]. The strong interaction of quercetin with DNA is responsible in part for its strong biological activity (in Chinese hamster V79 cells) since its dihydro counterpart (dihydroquercetin) with non-planar conformation did not bind to DNA and showed a weaker activity [9,10]. However, recent studies suggest that quercetin shows carcinogenic and mutagenic properties probably in part due to its pro-oxidant effects.
In this study we studied the interaction of dendritic antioxidants, a previously in-house made potent G0 dendrimer derived from syringaldehyde, G0 dendrimers derived from vanillin and 5-iodovanillin, with pBR322 DNA and compared with those of quercetin and vitamins, C and E using computer simulation.
The molecules with strong interaction with biomolecules may offer beneficial properties for many biomedical applications. Materials desired for biomedical applications require well-defined size, shape, and physico-chemical properties. Dendritic nanomolecules have unique architectures, which is beneficial for biological applications. In general, dendrimers are regular tree-like macromolecules accessible by chemical synthesis from a variety of building blocks. Their unique topology enforces a globular shape that has interior void space and multiple copies of functional groups at the termini of the dendritic branches. Therefore, dendrimers have been exploited as a drug delivery vehicle via encapsulation in the void space as well as covalent attachment on the surface [11,12]. Multiple surface functional groups can have useful applications where one needs signal amplifications. Considering the advantageous aspects of dendrimer, we have designed antioxidant molecules in a dendritic form for potential biological applications. The dendritic antioxidants have a phenolic hydroxyl group at terminal of each branch and benzylic hydrogen within the branch. The phenolic hydroxyl groups and benzylic hydrogen will serve as free radical scavenging points.
Although experimental methodology and the results are far superior to computational methods, it involves the arduous experimental procedures, thus labor-intensive and expensive. With recent advances that have increased the speed of computations and with greater sophistication of software programs, it is now possible to undertake computational projects related to molecular interaction between drug and DNA or other biomolecules like protein and predict the outcomes. There are theoretical studies showing that DNA-drug intercalation can be well understood by simulation and have good agreement with experimental research [13,14]. Two computational methods, namely molecular mechanics modeling and molecular dynamics (MD) simulations, have provided significant insights into the structure properties. This research will make use of some of the theory and methodology related to these areas. Recent monographs onthis topic contain not only the physical and mathematical background necessary to use the method but also point to applications of the methods [15][16][17].
Molecular based understanding of biological system can be well understood by the atomistic based simulation which complements experiment. Simulation results at the atomic level interactions between molecules help to analyze and interpret the experimental data. Over the the years various simulation methods have been applied to many problems in biological systems. The results have shown enough evidence in analyzing how enzymes behave, and how biomolecules adopt their functional structures depending upon their operating environment. Simulation results can help in the drug design process as well as understanding the molecular basis of disease in biological system. Computer simulations tremendously helped us to understand the way biological molecules move, stretch, and contract under dynamics conditions. These are essential to their functioning and morphological appearances, which justifies Feynman's assertion sixty years ago.
Molecular recognition (MR) technique is often used to study the non-covalent interactions among two or more compounds belonging to host-guest family such as enzyme-and drug-receptor complexes [18,19]. Drug design validation for treatment and prevention is challenging and is ideal for computational method in order to understand the thermodynamic predictions of the molecular identification process. Herein, we report computational molecular dynamics techniques to study various complexes with an unbiased fashion. There are literatures available regarding the computational and experimental studies of DNA with various drug interactions [20][21][22]. With respect to recently reported computational studies on Cyclodextrin inclusion recognition mostly by molecular dynamics technique, which is based on a systematic, automatic and quasi-flexible docking approach that prevents the influence of the chemist's intuition in the configuration generation [23,24], the inclusion molecular system focused in this work is based on pBR322 as a host ( Figure 1) and a set of dendritic molecules and commonly used antioxidants as guests. In an effort to outline a chart of interaction of DNA-antioxidant energetic, we present here a detailed thermodynamic view of DNA (pBR322) minor-major groove recognition by the G0 dendritic molecules in comparison to commonly used antioxidants (quercetin, vitamins C and E).

Computational and simulation experimental methods
Molecular mechanics and molecular dynamics calculations were performed with the current state-of-the-art Discover and Forcite program [25] using the universal, consistent valence force field (CVFF) /dreiding force fields [26][27][28] on a PC workstation (HP, USA). Before the dynamic simulation, the energy of the dendritic antioxidant molecules, derived from syringaldehyde (SYR), vanillin (VNL), and 5-iodovanillin (IDO), and naturally available small molecule antioxidants, quercetin, vitamins C and E were minimized and the conformational search of antioxidants were performed by simulated annealing molecular dynamics-full energy minimization strategy [29]. The lowest energy conformations were obtained using the above method and production run for the simulation experiment was carried out again at least 2000 ps. Thestructures of these molecules were obtained by energy minimization of a crystallographic geometry [30]. Each of the antioxidants were interacted with pBR322 plasmid DNA with base pairs of 20 and 80 (chosen for this calculation and built based on the previously reported amino acid sequence) as a host with/ without solvents under various ionic conditions. All simulation data were obtained in water environment with specific pH condition and 1:1 stoichiometry. Analyses were conducted on energy minimized structures of pBR322, each antioxidant, and pBR322-antioxidant complexes as well as molecular dynamics trajectories developed for this study. Choice of the conformation of each ligand-pBR 322 DNA dependant complex was carried out by choosing the preferred complexes with free energy (ΔGc) lower than 1.0 kcal/mol, which is experimentally determined [31,32]. This level of binding energy is comparable to drug-receptor interaction complexes [31,32]. All the calculation was performed on HP Compaq PC 700 MHz.

Molecular mechanics and molecular dynamics principles
Molecular modeling of drug interactions in the active site of the protein/enzyme/DNA can offer helpful information to pre-estimate which drug molecules would be effective for a particular ensemble of resistance mutations by optimizing a treatment regime. The current paper focuses on antioxidants interaction and binding preferences on pBR322 DNA. Molecular dynamics simulations [33,34], indicate that average energy difference between the host-guest complexes, frequently used as a measure of chiral recognition to drug screening process, which mostly depends on the length of the simulation time. Recent advances have been focused on the refinement stages of the homology modeling such as the conformational prediction of DNA and surface side chains as well as the modeling of ligand/receptor induced fit effects [35,36] which are essential steps to make the model useful as a drug discovery and optimization target. These kinds of high-resolution hostguest structure prediction applications have generally been performed using atomistic physics-based free energy estimators using molecular mechanics and molecular dynamics simulations. Molecular mechanics is a mathematical formulation to compute molecular properties of large number of atoms by mimicking the molecular geometries, energies and other features by adjusting bond lengths, torsion and angles, to equilibrium conditions, but different from ab initio calculation. Such study will enable us to understand the effect of geometry, stability of conformers, various spectra and rotational barriers. By employing various numerical schemes we can optimize the system/s under investigation as a function of energy in order to obtain the a minimum energy and time-dependent behavior of molecules.

Monte carlo docking minimization simulations
The host and guest molecules were positioned in the neighborhood with a distance of 12-15 Å. Monte Carlo docking simulations started by conjugate gradient energy minimization of this initial configuration for 200 iterations and accepted it as the first frame. In the course of trial to a new configuration, guest molecules could take translational movement of maximum 3 Å to x, y and z axis and rotation of maximum 180° around x, y and z axis. For docking with guest molecule flexibly, the torsional angles of three single bonds could rotate up to 180°. There are 3translational, rotational and dihedral rotational movements each, which makes nine degrees of freedom for the present system. Each cycle began with a random change of at the most 5 degrees of freedom among them. If the energy of the resulting configuration was within 1200 kcal/mol from the last accepted one, it was subjected to the 200 iterations of conjugated gradient energy minimization. The energy tolerance of 1200 kcal/mol was imposed to avoid significant overlap of van der Waals radii in the random search. After the energy minimization, the resulting structure was accepted based on criteria of a root-mean-squared displacement (RMSD) check, which compared the RMSD of the new configuration against those accepted so far. Configurations within 0.2 Å RMSD of pre-existing ones were discarded to avoid accepting similar configurations. The Monte Carlo simulations were performed until energy convergence. A cut-off distance of 5 Å was imposed on the calculation of non-bonded interactions, and the dielectric constant was set to 4, which imitate a solvent environment and Boltzmann averages of energies were evaluated at 298°K. The lowest energy structures obtained from the MC docking simulations were used as starting conformations in all of our further molecular dynamics simulations, without imposing any cut-off distance for non-bonded interactions. Constant NVT molecular dynamics calculations were performed using the leap-frog algorithm with a 1 fs time step. The initial atomic velocities were assigned from a Gaussian distribution corresponding to a temperature of 298°K. The system was equilibrated for 1000 ps and production run was done for 4 ns. The temperature was controlled by velocity scaling in equilibration phase and by Berendsen algorithm [37] in production phases with a coupling constant of 0.2 ps. Intermediate structures were saved every 5 ps for analysis. The interaction energy was defined as nonbonding energies between guest and host molecules.

Results and Discussion
Packing and shape-size complementarities, Van der Waals and hydrophobic interactions of the dendritic antioxidant molecules in the minor-major groove of the pBR322 DNA in conformity studies revealed many interesting and favorable electrostatics by desolvation. Hydrogen bonds appeared to be important for specificity based on structural considerations. This is not always important in understanding the binding interactions on post analyses of molecular dynamics trajectories. It may be due to the thermal averaging, solvent, and counter-ion effects. The strength of the hydrogen bonds retained between the pBR322 and the dendritic ligands of syringaldehyde during the molecular dynamics simulations were approximately 1.1 and 1.4 kcal/moles, respectively. In the present study, the three above-mentioned molecularly diverse compounds were selected on the basis of the standard free energy of complexation with pBR322 and have been used as a guidance set of molecular docking and interaction for computer simulation experiments. Statistical and graphical analyses of the inclusion geometries and binding energy of complexation were the main focus of the simulation. In order to study interaction modes of antioxidants with DNA, we used pBR322, a circular plasmid DNA, as a model DNA. Since free radical can attack any segment of DNA and cause oxidative damage, we built a 20 base pairlong DNA based on the published nucleotide sequence. The DNA was further extended to build 80 base pair long DNA. Final energy minimized conformations of (a) 20 base pair and (b) 80 base pair pBR322 are shown in Figure 1.
As shown in Figure 2, the energy minimized conformations of dendritic antioxidants derived from syringaldehyde (a), vanillin (b), and 5-iodovanillin (c), and commonly used small antioxidant molecules, quercetin (d), vitamin C (e), and vitamin E (f) were determined. All of the lowest conformation of dendritic antioxidants had three dimensional conformations. In comparison, vitamins C and E had a planar structure with a hydrocarbon tail. Quercetin has a non-planar structure consisting of two rings fused together to be co-planar and catechol ring with torsional angle of 153°. Russo et al. [38] determined the structure of quercetin with AM1 calculations. In their study the lowest energy minimum conformation occurred at 153° and then next lowest minima at 27°, which are comparable to our results. The authors also stated that quercetin requires only 0.5 kcal/mol to acquire planar conformation, indicating that the energy difference between nonplanar (153°) and planar (180°) conformation is negligible.

Monte carlo docking minimization simulations
Modeling the intermolecular interactions between a molecule and a biomolecule e.g., drug and its receptor or enzyme, is complex problem. During the modeling processes one need to take into account of various intermolecular forces such as covalent and non-covalent bonding( dispersion, van der Waals, hydrogen bonding, hydrophobic and electrostatic interaction) [39,40]. Binding forces is due to hydrophobic interactions, but the specificity could be controlled by hydrogen bonding and electrostatic interactions. In addition to these the effect of solvent on the binding association and the degrees of freedom need to be taken into consideration. Therefore, Monte Carlo docking technique tries to mimic the natural course of interaction between the molecules until the lowest energy acquires. In the process of Monte Carlo (MC) docking simulation, antioxidants were considered as a guest-ligand and pBR322 as a host-receptor molecule. During the simulations, the whole coordinates of pBR322s were flexible. The pathways of MC docking simulations showed a general tendency of forming inclusion complex and lowering interaction energy. The interaction energy was defined as the difference between the sum of the energy of individual guest and host molecule and the energy of the inclusion complex [41,42]. We tried several MC runs to search for the lowest energy configuration of the complex. Figure 3 shows a typical low

Molecular dynamics simulations
One of the main objectives of this research is to specifically understand how the energy use and movement of plasmid-DNA is impacted by various dendritic antioxidants and other antioxidants which are used routinely in biomedical applications. Such information is vital to the understanding of currently utilized cancer drugs, with or without supplements such as antioxidants for breast, ovarian and lungs [43]. Therefore, we used the lowest energy structures from the MC docking simulations as starting conformations for further molecular dynamics simulations [1,44]. Figure 4 shows the final energy minimized structures of 80 base pair pBR322 DNA complexed with 1 SYR. The energy of the complex was 2.52e4 kcal/mol. Typical final configurations of 20 base pair pBR322 complexed with SYR, VNL, IDO, vitamin C, quercetin, and vitamin E are shown in Figures 5-10. In all cases, the antioxidant molecules moved towards the minor groove of the pBR322 DNA.
In the molecular dynamics simulation using three SYR molecule and 20-mer DNA, two SYR molecules sit in the minor grove of the 20mer DNA and they are interacting with phosphate groups via hydrogen bonding. However, one molecule was intercalated between the bases (Figure 5a). The bulky molecule disconnects Watson-Crick hydrogen bonding and forms H-bonds and hydrophobic interaction with the bases. In case of five SYR molecules, four out of five SYR molecules are associated with DNA through external binding (Figure 5b). The fifth molecule did not have space to occupy and simply hanging off the phosphate group. Binding energy of each complex was 2.90e3 kcal/mol and 6.64e3 kcal/mol for 3 SYR-20mer DNA and 5 SYR-20mer DNA, respectively.    Three VAN molecules are also externally associating with DNA, two VAN molecules sit alongside of DNA and one VAN make a contact with DNA through two periphery phenol rings (Figure 6a). All three molecules are nearby phosphate groups, indicating hydrogen bonding is a major interaction force. In the dynamic simulation using five VAN molecules, some of them are intercalating and the rest binds externally (Figure 6b). The energy of the each complex are 2.88e3 kcal for 3 VAN and 1.04e4 kcal/mol for 5 VAN molecules.
In three IDO molecule dynamic simulations, two IDO molecules intercalate between bases (Figure 7a). Two periphery benzene rings containing IDO groups are facing the bases, making π-π stacking as well as H-bonds. Due to its steric bulkiness, bases are pushed away, thus disconnecting hydrogen bonds between bases. Five IDO molecules overpopulate the 20-mer DNA, leaving some of the molecules unbound (Figure 7b). The minimum energy of each complex is 4.55e3 kcal/mol for 3 IDO and 6.49e3 kcal/mol for 5 IDO.
Hydroxyl groups on the tail part of vitamin C molecule make H-bonds with phosphate groups (Figure 8a). The enolic ring intercalates and its hydroxyl groups form H-bonds with bases (amino and carbonyl groups) and phosphate groups. Since vitamin C is a small molecule, we also did simulation with 15 molecules (Figure 8b). It shows that vitamin C molecules can intercalate fairly well into DNA and has less specificity towards DNA binding site; Vitamin C molecules bind to minor groove as well as major groove. The minimum energy of each complex is 2.71e3 kcal/mol for 3 vitamin C and 13.22e3 kcal/mol for 15 vitamin C.
The major interaction force of quercetin to DNA determined using three (Figure 9a) and five quercetin molecules was hydrogen bonding formed with phosphate groups. This study shows that four out of five quercetin molecules bind externally and only one molecule intercalates (Figure 9b). Figure 9 also shows that before (c) and after (d) energy minimization of 20mer-pBR322 complexed with ten quercetin molecules in the dynamic simulation experiment. Due to semi-planar structure of quercetin, it was expected to intercalate fairly well. However, surprisingly, only a few quercetin molecules intercalate through the catechol ring and most of them are bound externally, forming H-bonds. The minimum energy of each complex is 2.76e3 kcal/ mol, 3.89 e3kcal/mol, and 32.18 e3kcal/mol for 3, 5, and 10 quercetin molecules, respectively. Surprisingly, vitamin E binds well with DNA and intercalates through phenolic moiety, forming π-π stacking interaction ( Figure  10b). Figure 10 also shows the energy minimized pBR322 associated with three (c) and five (d) vitamin E molecules. Most of vitamin E molecules bind DNA through externally but some of them are able to destroy Watson-crick base pair and confine themselves between the two separated chains. This experiment suggests that hydrophobic vitamin E can bind DNA well although it is questionable as to how effectively it can arrive nucleus in aqueous environment. Figure 11 shows the binding energy profiles of 20-mer pBR322 complexed with various number of antioxidant molecules. The binding energy indicates their capacity to form thermodynamically stable inclusion complex with DNA. Therefore, from the energy profile of each antioxidant, we can determine the effective number of the antioxidant required to attain the minimum energy complexation with DNA for a given number of base pair. Based on the graphs, SYR forms thermodynamically more stable complex with DNA than any of the antioxidant molecules we used in our computer simulation experiment. It appears that its lowest binding energy was achieved at -70 kcal/mol with 3 SYR molecules for 20 base pair long DNA and at -63 kcal/mol with 10 SYR for the same DNA. In other words, energy minima of SYR occurs somewhere between 3 to 10 SYR molecules for 20-mer pBR322 DNA. In comparison, quercetin, which showed slightly higher binding energy than SYR, had its lowest binding energy at -58 kcal/mol with 9 molecules. In order to see whether the number of base pair unit of pBR322 DNA has any effect on the total number of antioxidant molecules required to lower the binding energy, we also performed the same experiment with 80 base pair long-DNA (data not shown except SYR, Figure 12). Figure 12b shows that the binding energy of SYR-80-mer pBR322 reached its lowest optimum at -800 kcal/mol with around 15 to 20 SYR molecules, implying that as the base pair of the DNA increases to 80, the ratio between the number of antioxidant and base pair increases to approximately 0.25:1. Due to computer memory constraints, we could not carry out further studies with increasing number of base pairs. It will be continued in the future using high-power cluster network.
In terms of binding energy of antioxidant to pBR322 DNA, SYR formed the more stable inclusion complex with pBR322 irrespective of the total number of bases; SYR can achieve the lower binding energy with the fewer number than other antioxidants. If the number of antioxidants associated with DNA at lowest energy is extrapolated to overall efficiency of the antioxidant, SYR is the most efficient antioxidant among the antioxidants used in this simulation experiment and that SYR might be more effective and promising antioxidant than quercetin [45] and vitamin E which are being considered to be effective against a wide variety of inflammatory diseases, including cancer. Although quercetin is known to produce various beneficial biological effects [46,47], many recent studies have reported that repeated usage of quercetin can be also carcinogenic and mutagenic, especially in the presence of transition metal ions like copper ion. Some studies are currently undergoing regarding its safety and efficacy against various diseases.
Although all antioxidant molecules used in this study moved towards the minor groove of the pBR322 DNA, some moved faster than others. Therefore, we determined mean square displacement (MSD) of SYR and quercetin with respect to DNA since they had the lowest binding energy among the antioxidants we used in this study. The MSD data provide useful information as to which site of DNA antioxidants bind to as well as how fast the antioxidants find their favored site. Three molecules of each antioxidant were used for this experiment. It showed that both quercetin and SYR prefer to be near the guanine-cytosine close to the phosphate group of the minor groove as shown in Figure 13. Typical MSDs for SYR and quercetin are shown in Figure 14 and Figure 15, respectively. These graphs illustrate how fast these antioxidants move towards the minor groove of the pBR322 DNA under the right condition and structural interaction of the host to the guest. SYR molecules reached to the minor groove less than 49 ps and remained near the site. In the case of quercetin, one of the three quercetin molecules (top curve in Figure 15) reached inside of the minor groove of the 20-mer pBR322 DNA by 28 ps and its mean square displacement is around 241 Å2. However, within a short period of time, the quercetin molecule moved away from the DNA site, showing extremely unstable behavior during the study of our simulation experimental time. The second and third quercetin molecules showed  oscillatory motion which is also indicative of instability of the quercetin-DNA complex and such motion was not observed with SYR. The mean square displacement data of quercetin probably mean quercetin might be a short acting antioxidant.
Taken together, it is clear that dendritic SYR antioxidant molecule is the more preferred structure for the host pBR322 than quercetin computationally. The ability of SYR to interact well with DNA might make it potent antioxidant to guard DNA from harmful oxidative damages, eventually lowering the risk of mutagenesis and carcinogenesis, which of course need further investigation theoretically and experimentally. Based on this study, we conjecture that SYR might better protect DNA than quercetin.
Some of the antioxidants like vitamins C and E have a chiral center. In this study, their chiral interaction was ignored. The prediction of chiral recognition at a molecular level will be challenging but valuable in comparing binding modes of such antioxidants with non-chiral antioxidants and evaluating their antioxidant activities. Generalized molecular modeling methods for enantio-discrimination by other various dendritic antioxidants are under investigation.
Since the effects of proteins, salts, and divalent cations which coat cellular DNA are still unknown, the antioxidant molecules were bound to a naked pBR322 DNA plasmid in this study. It should therefore be noted that this system is still far from an in vivo model and the calculated results need to be further validated under various physiological conditions, by different characterization techniques if and when available, and possibly by wet-lab experiments.

Conclusion
Molecular modeling studies have given us useful insight into interaction of various antioxidants with the plasmid DNA at the atomic level. Throughout this research, MD simulation was carried out to investigate the effects of antioxidant complexation with pBR322 in terms of the relative distribution in the interaction and binding energies.
Our study suggests that all of the antioxidants including the dendritic antioxidants reside near the minor groove of guanine-cytosine site of pBR322 DNA. Interaction energy profile indicated that the antioxidant SYR formed more favorable interaction with pBR322 DNA than all other antioxidants used under this study. To our knowledge, this is the first report of its kind on antioxidant-binding specificities on a native DNA molecule using such low concentrations. As a first hand study we tried to enhance our understanding on interaction and binding energy of antioxidant with pBR322 DNA by doing computer simulation and established computational strategies in order for experimentalist to determine if further investigation will be pursued in this case.