Comparative antibacterial analysis of the anthraquinone compounds based on the AIM theory, molecular docking, and dynamics simulation analysis

Hydroxyanthraquinones and anthraquinone glucoside derivatives are always considered as the active antibacterial components. Comparison of structure characteristics and antibacterial effect of these compounds was performed by applying quantum chemical calculations, atoms in molecules theory, molecular docking, and dynamics simulation procedure. Density functional theory calculation with B3LYP using 6-31G (d, p) basis set has been used to determine ground state molecular geometries. The molecular geometric stability, electrostatic potential, frontier orbital energies, and topological properties were analyzed at the active site. Once glucose ring is introduced into the hydroxyanthraquinone rings, almost all of the positive molecular potentials are distributed among the hydroxyl hydrogen atoms of the glucose rings. In addition, low electron density ρ (r) and positive Laplacian value of the O–H bond of the anthraquinone glucoside are the evidences of the highly polarized and covalently decreased bonding interactions. The anthraquinone glucoside compounds have generally higher intermolecular binding energies than the corresponding aglycones due to the strong interaction between the glucose rings and the surrounding amino acids. Molecular dynamics simulations further explored the stability and dynamic behavior of the anthraquinone compound and protein complexes through RMSD, RMSF, SASA, and Rg. The type of carboxyl, hydroxyl, and hydroxymethyl groups on phenyl ring and the substituent glucose rings is important to the interactions with the topoisomerase type II enzyme DNA gyrase B.


Introduction
With the widespread use of antibiotics and the emergence of resistant bacteria, there is a need for developing antibacterial compounds that have new mechanism of actions against bacterial strains. Therefore, the search for new antibacterial agents remains a significant challenge in the future. Natural products are structurally complex and pharmacologically valuable products; the separation, synthesis, and the structure-activity relationship attract extensive attention for scientific research. There are three authorized Rhubarbs in Chinese Pharmacopoeia, namely, Rheum palmatum L., Rheum tanguticum Maxim. ex Balf, and Rheum officinale Bail, with the roots and rhizomes generally called official rhubarbs. Rhubarb contains a lot of active compounds, such as anthraquinone, anthrone, butyrophenone, tannin, and some glycoside derivatives, and it also has the functions of purging heat, clearing the intestines, cooling blood, lowering blood sugar, antibacterial, antiviral, antitumor, and lipidlowering [1]. The hydroxyanthraquinones are considered as the active antibacterial components, such as aloe-emodin, rhein, emodin, and chrysophanol [2]. In addition, some anthraquinone glucoside derivatives are also attracting more and more attentions due to their good microbicidal activities [3,4]. It was reported that the anthraquinone glucosides were much more potent than their corresponding aglycones. For example, emodin has lower bacterial neuraminidase inhibitory activity than the corresponding glucoside derivative [5]. However, there are few comparative researches about the anthraquinone and its glycosides.
DNA gyrase is a member of topoisomerase type II enzyme, which controls the topology of DNA in the time of transcription, replication, and recombination processes [6]. DNA gyrase is composed of two A subunits (GyrA) and two B subunits (GyrB), where the GyrA is involved in DNA transit and the GyrB subunit contains an ATP-binding domain [7]. DNA gyrase has proven to be a worthy target for antimicrobial agents [8]. A lot of new inhibitors were developed and evaluated against the DNA gyrase B to explore new antibacterial drugs [9][10][11]. Due to the absence of detailed report about the relationship of the structure-antibacterial activity for these aglycones and their anthraquinone glucoside derivatives by now, the aim of present text is to give a comprehensive comparison of these compounds derived from rhubarb by applying standard quantum chemical calculations [12], Bader's theory [13] of atoms in molecules (AIM), molecular docking procedure [14][15][16], and dynamics analysis selected the DNA gyrase as a target.

Computational methods
The molecular structure optimization and corresponding energies of the compounds are calculated by using GAUSS-IAN 09 W program package [17]. It was found that density functional theories (DFT) with Becke-3-Yang-Parr (B3LYP) lead to the most successful correlations between observed and calculated vibrational wavenumbers and have emerged as a good compromise between computational cost, coverage, and accuracy of results [18,19]. The DFT with B3LYP combined with standard basis sets 6-31G (d, p) was thus applied to optimize these anthraquinone compounds in this text. Topology analysis of electron density and the molecular electrostatic potential surfaces on van der Waals (vdW) surface in both forms were performed by using the Multiwfn program [20]. The protein of the template structure DNA gyrase B (PDBID: 4DUH) was downloaded from PDB website (http:// www. rcsb. org). In target DNA gyrase B preparation, the protein was prepared by using the biopolymer module, through which the missing atoms were identified and the repair sidechain was repaired. The crystal water (H 2 O) molecules and other ligands were deleted from the crystal structure. The polar hydrogen atoms and Gasteiger charge were added. The protein molecule was also assigned AD4 type. The ligand was docked into the catalytic pocket of the DNA gyrase B based on the AutoDock 4.2 program [21]. 4-{[4′-Methyl-2′-(propanoylamino)-4,5′-bi-1,3-thiazol-2-yl]amino}benzoic acid was considered as the reference ligand. The Gasteiger charges of the receptor protein and the anthraquinone compounds were computed and saved as PDBQT format. The genetic algorithm with default settings was selected, and the number of runs was set to 100. The configuration file with binding site information was generated with grid coordinates 3.024, 2.008, and 36.865 for x -, y -, and z-axis, respectively. The docking grid box was set at 30 × 30 × 30 points with a spacing value of 0.375 Å. Genetic algorithm with default settings was employed, and the number of runs was set at 100. Intermolecular interactions exist between these anthraquinone compounds, and the amino acids residues, such as Asn46, Asp73, Ile78, Thr165, Arg76, Ile9, Gly101, Lys103, Val120, and Arg136, play an important role in their stability at the binding site. Topological analysis was performed by using AIMALL package [22], and the corresponding wave functions were also generated in order to gain further insights into the bonding characters.
In order to evaluate the stability of obtained compounds at the active site, molecular dynamics (MD) simulations were carried out by using Gromacs 5.1.4 software. The ligands were simulated parametrically with the AMBER force field GAFF by using the online program [23]. The protein topology and coordinate files were generated through the Amber 99SB force field [24]. The solvent box of this system was a regular dodecahedron filled with SPC/E water molecules [25], and the distance between the compound and the edge of the box was 10 Å. Na + ions were then added to ensure the electrical neutrality system. The energy minimization of the complex was achieved by the steepest descent minimization algorithm. Then, the system performed isovolumetric isothermal (NVT) ensemble equilibration for 500 ps and used a 2 fs time step, followed by isothermal isobaric (NPT) dynamics at Langevin thermostat (310 K) and Langevin barostat (1 atm.), respectively [26]. Long-range electrostatics was calculated by using the particle grid Ewald method, and the temperature and pressure coupling was performed through a modified Berendsen thermostat and Parrinello-Rahman barostat. Analysis of the result was carried out by using Gromacs tools [27].

Molecular geometry and docking
According to the literature [28,29], the binding modes and binding affinities play important roles in the inhibition of 1 3 the small molecules to the enzyme. These anthraquinone compounds were firstly subjected to geometry optimization at the ground state; further molecular docking results indicated that the selected binding affinities for these complexes of molecule-DNA gyrase B are from − 7.27 to − 10.24 kcal/ mol. It was worth noting that the bonding energies of the anthraquinones are − 8.50, − 8.24, − 7.27, − 7.91, and − 7.66 kcal/ mol for molecules M1-M5, respectively. However, the corresponding anthraquinone glucoside compounds have generally higher intermolecular binding energy due to the strong interactions between the glucose ring and the surrounding amino acids Asn46, Arg76, Pro79, Gly101, and Lys103 (− 1 0.24, − 9.47, − 9.74, − 9.31, and − 9.11 kcal/mol for molecules M6-M10, respectively). These bonding energies are lower than those of ciprofloxacin (− 7.96 kcal/mol), which is used as a positive control. This is consistent with the result that the anthraquinone glucosides have much more antibacterial effect than their corresponding aglycones [30]. The optimized conformation of the most active compounds M1 (emodin) and M6 (emodin-8-beta-glucoside) with atom numbering scheme and the optimal docked postures and molecular interactions at the active site of DNA gyrase B are shown in Fig. 1. The residues Asn46, Asp73, Gly77, Ile94, Gly101, and Arg136 played an important role in the interactions in the 4DUH active site [31]. Comparative analysis showed that geometric conformations of these compounds changed obviously after interacting with the DNA gyrase B at the active site (Fig. S1).
There are obvious changes about the bond angles and dihedral angles rather than the bond lengths, which suggested strong intermolecular interactions. It can be seen from Table S1 that there is almost no change in the length of each compound bond at the active site. The glucose ring that binds to the aromatic ring ketone is obviously bent at the active site, which can be ascribed to the more serious intermolecular interaction with surrounding amino acids [23], such as Asn46, Gly101, Ala100, Ile94, and Phe104.
Compared with the compound M6, the glucose ring of M7 at the active site rotates clearly clockwise, resulting in a smaller angle with the anthraquinone ring plane. The dihedral angle of C11-C12-O18-C26 is changed from 16.23° at the optimized state to − 116.95°, which can be ascribed to the strong hydrogen bonding and electrostatic interactions with surrounding amino acids Asn46, Ala100, Phe104, Gly117, and Val120. Besides the multiple hydroxyl groups on the anthraquinone skeleton, hydroxyl groups of the glucose ring are also important to the antimicrobial effect, which may be ascribed as the hydrogen bonding interactions at the active site. In fact, the combined anthraquinone chrysophanol 8-O-beta-D-glucoside shows good bioactivity compared with the free chrysophanol [32]. Therefore, the influence of substituent groups on anti-DNA gyrase B activity might be related with not only the type and number of substituent carboxyl, hydroxyl, and hydroxylmethyl on phenyl ring groups but also the combined glucose rings. The nearest neighbors and short contact distances (Å) of anthraquinone compounds with the amino acid residues of DNA gyrase B are shown in Table 1.

Molecular electrostatic potential
Molecular electrostatic potential (MEP) on molecular vdW surface is an important parameter to study the intermolecular interaction [33][34][35]. Usually, MEP is also related to the electronic density [36]; thus, it is always used to anticipate reactive sites that electrophilic and nucleophilic reagents attack [37]. In this text, the molecular electrostatic potential of the compounds at the active site was compared to investigate the interactions between the complexes. In Fig. 2, the negative (green) and positive (red) regions are regarded as electrophilic and nucleophilic reactivity regions, respectively, while the white regions show zero potential. In general, the electrostatic potential distribution of anthraquinone compounds is more uniform than that of the anthraquinone glucosides. The main positive charge is distributed near the hydrogen atom of the hydroxyl group on the side chain of the aromatic ring plane. For the molecule M3, the oxygen atom O18 in the carboxyl group acts as an electron donor to form a hydrogen bonding interaction with the HH11 of the amino acid Arg136, and thence, the region at the vicinity of the connected hydrogen atom H28 remains the highest electrostatic potential value of 66.73 kcal/mol. Due to the strong interaction between the side chain hydroxyl oxygen atoms (O17 and O19) of molecule M1 and the surrounding amino acid Gly77 and the electrostatic interaction between the oxygen atom O15 of the anthraquinone ring and amino acid Arg76, the electrostatic potential has reached the lower value of − 58.25 kcal/ mol (Fig. 2a).
Once glucose ring is introduced into the hydroxyanthraquinone plane, it was shown that almost all of the positive molecular potentials of anthraquinone glucosides are distributed among the hydroxyl hydrogen atoms of the glucose structures except the molecule M9. These substituent hydroxyl groups on the molecular aromatic ring interact with the amino acids Glu42, Asn46, and Gly11 by hydrogen bonding and electrostatic interaction, which might lead to imbalance of distribution of electrons (the value of the global surface maximum is 75.40 kcal/mol, Fig. 2i). In general, due to the presence of the hydroxyl group in the glucose ring, it can interact with the surrounding amino acids Lys103, Asn46, Ala100, Phe104, and Val120 by intermolecular hydrogen bonds, electrostatic attraction, and van der Waals forces, resulting in major changes in the molecular structure and the molecular electrostatic potential distribution. In general, the positive molecular potentials at the vicinity of the molecular hydrogen atom are increased obviously while glucose ring is introduced into the hydroxyanthraquinone ring. For example, the global surface maximum of M1 is 61.35 kcal/mol while the value changes to be 70.18 kcal/mol for M6, and the value of M2 (44.37 kcal/ mol) increases about 8.19 kcal/mol (52.56 kcal/mol) for M8. This indicated that the electronic cloud of the aromatic side chain has shifted while interacting with the bonding site.

The molecular orbitals
The important frontier molecular orbitals (FMOs) are the highest (HOMO) occupied molecular orbitals and the lowest (LUMO) unoccupied one, which play an important part in the chemical stability of the molecule [38]. The HOMO represents the ability to donate an electron, while the electron acceptor LUMO represents the ability to accept an electron. The energy gap between HOMO and LUMO also determines the chemical reactivity, optical polarizability, and chemical hardness-softness of the molecule [39]. According to Koopmans' theorem, E HOMO and E LUMO values of any chemical type are associated with its ionization energy and electron affinity values, respectively [40]. The energy differences between the LUMO and HOMO are known as energy band gap ΔE . In addition, for closed-shell compounds, hardness ( ), chemical potential ( ), and electronegativity ( ) are also defined. Recently, a new descriptor electrophilicity index ( ) has also been defined to quantify the global electrophilic power of the compound [41]. Dipole moment D (Debye) is also shown in Table 2. The usefulness of this new reactivity quantity has been recently demonstrated in understanding the toxicity of various pollutants in terms of their reactivity and site selectivity [42].
where A and I are the ionization potential and electron affinity of the compounds, respectively.
Commonly, the atom occupied by more densities of HOMO should have stronger ability to detach an electron whereas the atom with more occupation of LUMO should have ability to gain an electron. In the present study, the HOMO and LUMO energies are predicted at B3LYP method with 6-31G (d, p) basis sets. It is clear from Fig. 3 that the positive and negative phase is represented in orange and green color, respectively. It can be seen from the plots that the region of HOMO and LUMO levels spreads over the free anthraquinones (M1-M5). However, the HOMO and the LUMO orbitals are located all over the anthraquinone ring plane of anthraquinone glucosides rather than the glucose ring. This indicated that the atom with more occupation of LUMO should have ability to gain an electron [43]. However, the hydroxyl oxygen atoms in the glucose ring act as electron donors and easily form hydrogen bonding interactions with the surrounding amino acids, resulting in stronger ability to detach electrons. It is interesting to find that the glucose ring donates electrons for M9, while the anthraquinone ring plane accepts electrons. From the energy levels of HOMO and LUMO orbital for these anthraquinone compounds at the active site (Table 2), it can be seen that the energy gap of HOMO-LUMO is the highest for M9 (0.148 eV), which reflects the chemical stability of the molecule.
Due to the fact that molecular resistance can change or deform the number of electrons, chemical hardness ( ) is associated with the stability and reactivity of a chemical system [44]. The combined anthraquinone M9 has the highest ∆E gap , which suggests a stable molecule and it needs much excitation energy to reach the manifolds of excited states [39]. As can be seen, the calculation electrophilicity ( ) performed at anthraquinone glucosides is very close to each other. The other important electronic property dipole moment values for the molecules are also given in Table 2. Generally, the combined anthraquinone glucosides have larger dipole moments than the free anthraquinones.

Topological property analysis
In general, the charge concentration/depletion at the BCP (bond critical point) of the chemical bond provides interesting information about the chemical bonds by using the Laplacian of electron density ∇ 2 ρ bcp [45,46]. The electron density (r) and its Laplacian ∇ 2 (r) help to determine the nature of interactions. In general, large values of the density (r) and its Laplacian ∇ 2 (r) show the power of bond interactions. Negative values of ∇ 2 (r) indicate a strong covalent character, whereas the positive one indicates a decreased charge in the internuclear region [47,48]. Due to the fact that emodin (M1) and its anthraquinone glucoside emodin-8-beta-glucoside (M6) have the  [18]. It can be seen from Fig. 4 that the electron densities of the molecule M6 at active site are very close to the optimized phase. However, the electron densities of the M1 at active site are higher than those of the other status. The electron densities (ρ bcp ) of the C-C bonds on aromatic ring for these molecules in both forms ranged from 1.583 to 2.159 eÅ −3 , implying strong covalent interactions [44]. Due to the strong electrostatic interactions between the carbonyl O7 atom and amino acids Asn46, His99, Phe104, and Val118, the electron density of the bond C6-O7 remains the highest value. In short, almost all of the electron densities of M1 are reduced at the active site, while the M6 remains constant (Fig. 4). Compared with other chemical bonds, the electron densities of all O-H and C-O bonds are generally higher while the corresponding Laplacian values are lower except C6-O7, which suggests that the charge densities at the BCP are highly concentrated.
On the other hand, low value of electron density and positive character of the Laplacian for M1 are the evidences for the covalent rupture of C6-O7 bond. Moreover, a change in the sign of the Laplacian from the positive values (4.952 eÅ −5 ) to a negative character at the active site (− 6.037 eÅ −5 ) indicates the depletion of the charge density (C6-O7 bond of M1) when it interacts with the DNA gyrase B. In general, most of the Laplacian values are less than 0, indicating the strong covalent interactions in the molecules (Fig. 4), which is different from the weak interaction of hydrogen bonds reported in the literature [47]. Once glucoside ring is introduced into the M1 (Fig. 4), low value of electron density (r) and positive data of the Laplacian of the O19-H27 are the evidences of the highly polarized and covalently decreased bonding interaction [49]. Actually, most of the Laplacian values of electron densities are reduced at the active site. Among which, the highest electron density values (2.766-2.788 eÅ −3 ) can be found on the carbonyl C-O bonds of the parent ring in both situations. Most of the electronic Laplace values for M1 molecule became more negative on the anthraquinone ring plane, suggesting a more concentrated character at the active site. The electron densities of the three hydroxyl groups on the M1 hydroxyanthraquinone plane become more depleted after interacting with the DNA gyrase B, implying strong hydrogen bonding interactions with the surrounding amino acids Gly77, Gly101, and Asp25 [50]. Most of the covalent bonds of the M6 molecule become more

Molecular dynamics (MD) simulation
In order to further explore the dynamic behavior of the complexes M1-DNA gyrase B (M1-GB) and M6-DNA gyrase B (M6-GB), MD simulation was performed. The stability of complex was evaluated by changes in root mean square deviation (RMSD) and root mean square fluctuation (RMSF) as well as changes in the system energy. As shown in Fig. 5a, the RMSD curves of the M1-GB and M6-GB complexes are similar to each other, and they reach stability around 20 ns. The final RMSD ranged from 1.3 to 1.7 Å, and the curves tend to be flat, indicating equilibrium systems [51]. The RMSF of the two complexes are similar, suggesting similar residues that play a key role in the interaction of the two complexes (Fig. 6b). The radius of gyration (Rg) is usually used to evaluate the compactness of system, and it indicates the volume and shape of the protein. As shown in Fig. 6c, the Rg value of M1 is slightly reduced at 0-5 ns, and it becomes stable after 5 ns. However, the Rg value of the M6-GB complex is increased slightly at the first 5 ns, and then, the Rg is gradually decreased, indicating that the system becomes more stable [52]. The solvent accessible surface area (SASA) is also calculated for the two studied complexes during 50 ns of molecular dynamics simulation. The SASA of the two complexes are found to slightly decrease, indicating that the anthraquinone compounds can affect the surface hydrophobicity, and the systems became more stable (Fig. 5d).  Figure 6a, b shows the number of hydrogen bonds between anthraquinone compounds and DNA gyrase B during the dynamics simulation. Analysis of the results indicates that the M1 could form 1 to 3 hydrogen bonds in average for whole 50 ns at the active site. On the other hand, M6 forms 2 to 7 hydrogen bonds in average at the active site of DNA gyrase B for the whole 50 ns of molecular dynamics simulation. This result confirmed the stronger molecular interaction between M6 and the protein, which was consistent with the molecular docking results. These results also suggested that the hydrogen bonds mainly participated in the interaction between the anthraquinone compounds and DNA gyrase B. In addition, the superpositions of docked conformations of the complexes before molecular dynamics and that after 50 ns of molecular dynamics are shown in Fig. 6c, d. It can be seen that the ligands retain the same binding sites before and after molecular dynamics simulation; the most evident conformational variation is the glucose ring, which plays a critical role in the molecular stability by molecular hydrogen bonding interactions with surrounding amino acids.

Conclusion
It was convinced that anthraquinone and anthraquinone glucosides were antibacterial active compounds [53]. It was concluded from this text that almost all of the anthraquinone skeleton planes are constant when interacting with the target protein DNA gyrase B. Compared with the multiple hydroxyl groups on the anthraquinone skeleton, hydroxyl groups of the glucose ring are also important to the antimicrobial effect, which may be ascribed as the hydrogen bonding interactions at the active site. It was consistent with those reported compounds in the literature that glycoside compounds showed moderate antimicrobial activities against Gram-positive bacteria and Gram-negative bacteria [54,55]. In contrast to the average distribution of anthraquinone compounds, MEP analysis indicated that almost all of the positive molecular potentials of anthraquinone glucosides are distributed among the hydroxyl hydrogen atoms of the glucose structures. In fact, the electron densities of hydroxyanthraquinone are generally higher than those of the corresponding anthraquinone glucosides. The electron densities become more depleted after docking with the DNA gyrase B, especially the hydroxyl groups. The molecular dynamics simulations further explored the stability and dynamic behavior of the complexes. The docked ligands retained the same binding sites before and after molecular dynamics simulation, suggesting the stable complex systems. It may be concluded that not only the functional carboxyl, hydroxyl, and hydroxylmethyl groups but also the glucose ring on phenyl ring can improve the antimicrobial activity. This indicates that glycosides played an important role in the antibacterial activity of natural products [56]. However, more experimental researches are needed to analyze the antibacterial structure-activity relationship, which will be useful to discover and design novel compounds.