Computational Binding Analysis of Ethyl 3,3,5,5-Tetracyano-2-Hydroxy-2-Methyl-4,6-Diphenylcyclohexane-1-Carboxylate in Calf Thymus DNA

In the present paper, several computational binding analyses were performed on ethyl 3,3,5,5-tetracyano-2-hydroxy-2-methyl-4,6-diphenylcyclohexane-1-carboxylate which was newly synthesized by three-component condensation of benzaldehyde with ethyl acetoacetate and malononitrile in the presence of trichloroacetic acid, and the structure was finally proved by X-ray analysis. The visualization of molecular interaction was carried out through Hirshfeld surface analysis and ESP. The atomic charges, HOMO, LUMO, and electrostatic potential were also studied to explore the insight of the molecule deeper, and then, natural bonding orbitals (NBO) and non-linear optical properties (NLO) were calculated to reveal the interactions that happen to be between the filled and vacant orbitals. Afterwards, molecular docking studies predicted the compound binding mode fits in the minor groove of DNA and remained interacts via stable bonding as validated by molecular dynamics simulations. The binding energy estimation also affirmed domination van der Waals and electrostatic energies. Lastly, the compound was found as good drug-like molecule and had good pharmacokinetic profile with exception of toxic moieties.


Introduction
Nowadays, the interaction studies of biomolecules are given importance to develop highly efficient target drug molecules. Binding sites are the places where molecular interactions occur. Hence, the analysis of binding sites is of crucial importance to understand the biological processes those various biomolecules are involved in. In order to derive the data responding to such scientific questions, we have carried out several steps, such as Hirshfeld surface analysis, natural bonding orbitals (NBO) and non-linear optical properties (NLO), and molecular docking to comprehend the molecular interactions and possible function of the selected compound. Such analysis has been performed on pyran-based heterocyclic compounds because they are having an irreplaceable place due to their distinct structures and great potential for further transformations. Pyran derivatives occupy an important position in the series of heterocyclic compounds of both natural and synthetic origins. Among those compounds, functionalized 2-amino-4H-pyrans containing electron-withdrawing groups in the 3-position constitute one of the most extensively studied subclasses. The presence of such functional groups as amino, cyano, or carbonyl in molecules of known 2-amino-4H-pyrans makes them promising intermediate products for the synthesis of fused heterocycles. Many fused 2-amino-4H-pyrans were found to exhibit biological activity. Compound shaving the pyran structural motif exhibits a wide range of biological activities, such as diuretic, analgesic, myorelaxant, anti-coagulant, anti-cancer, anti-tumor, and anti-HIV. In addition, they are also worth to be studied because of their benefits for the treatment of neurodegenerative disorders including Alzheimer's disease, amyotrophic lateral sclerosis, Huntington's disease, and Parkinson's disease [1][2][3][4][5].

Material
Herein, we perform several methods to discover the insight molecular interactions and possible binding modes and they have been performed on the three-component condensation of benzaldehyde with ethyl acetoacetate and malononitrile in the presence of trichloroacetic acid. It was found that, depending on the conditions, the condensation product was the expected product, which is so-called ethyl 3,3,5,5-tetracyano-2-hydroxy-2-methyl-4,6-diphenylcyclohexane-1-carboxylate (I). A single crystal suitable for X-ray analysis ( Fig. 1) was derived by double recrystallization from ethanol [6]. Crystallographic data and parameters of X-ray diffraction experiments are presented in Table 1. The complete set of X-ray diffraction data for compound (I) was deposited to the Cambridge Crystallographic Data Centre (CCDC entry no. 1 839 026) and is available at www. ccdc. cam. ac. uk.

Hirshfeld Surface and ESP Analysis
The Hirshfeld surface map helps to visualize the molecular interactions which take a significant role in the self-assembly process of a crystal. The conventional mapping of d norm and 2D fingerprint plots of newly synthesized three different heterocyclic compounds are drawn by Crystal Explorer 17.5 software [7] with the aid of Crystallographic Information File (CIF). And, the cluster of molecules within a 3.8 Å radius was selected around a reference molecule to calculate pairwise interaction energies in a crystal, and energy frameworks were also performed using CE-B3LYP/6-31G(d,p) technique.

Computational Details
The geometry optimization of the newly synthesized quinoxaline derivative has been carried out using Gaussian 09 at the B3LYP/6-311G** level of theory [8]. The corresponding threshold maximum force and displacement values are 0.00045 and 0.00018 au respectively. The atomic charges, HOMO, LUMO, and electrostatic potential map have been analyzed. The graphic illustration of iso-surface with a value equal to 0.2 of HOMO and LUMO was used. In order to understand the interactions that take a place among the filled and vacant orbitals, natural bonding orbitals (NBO) have been carried out. The graphical interpretation of all these studies was made by Gauss view and 3plot program [9,10].

Docking Study
Molecular docking of the compound with the DNA receptor (PDB id: 1BNA) was performed using AutoDock 4.2.6 software [11]. The DNA structure was retrieved from the protein databank and processed in UCSF Chimera [12], where co-crystalized water molecules were removed from the structure. The compound structure was considered for drawing the compound 2D structure, which was then converted into 3D coordinates and saved as.pdb format. Ligand preparation was done using Fig. 1 Structure of the molecule of ethyl 3,3,5,5-tetracyano-2-hydroxy-2-methyl-4,6-diphenylcyclohexane-1-carboxylate (I) according to the X-ray diffraction data AutoDock Tool (ADT), Gastegier charges were allocated, non-polar hydrogens were added, and then saved in PDBQ. The receptor DNA was subjected to the same processing as discussed above for the ligand molecule. Parametrization for the grid and docking calculations were conducted via ADT. The docking was performed by considering ligand as flexible while receptor as rigid. Blind docking was carried out by setting grid box size of 60 × 80 × 114 Å 3 with spacing of 0.375 Å. Eventually, the docking runs were accomplished using Lamarckian genetic algorithm (LGA) and an empirical-free energy function and at 100 binding conformation of the ligand molecule was evaluated. The lowest binding free energy binding mode the compound was selected and complexed with the receptor for post-docking analysis.

Molecular Dynamics Simulations and Binding Free Energies
The best-docked complex was further subjected to molecular dynamics simulation using Amber simulation package 20 [13]. The complex was submerged in a TIP3P water box, in which Na + counter ions were to neutralize the system. The OL3 force field was employed to generate parameters for the DNA molecule while GAFF [14] was used to parameterize the compound. The system was then energy minimized using combined action of the steepest descent and conjugate gradient algorithms. During the simulation process, NPT and NVT ensembles were applied to the system. SHAKE algorithm [15] was implemented to constrain bonds involving hydrogen atoms. The simulation production run was performed for 100 ns and the trajectories were processed through the CPP-TRAJ module [16]. The trajectories were then analyzed using the MMGBSA method [17]. In total, 500 frames were picked up from the simulation trajectories and subjected to the following MMGBSA equation for net binding free energy.

Drug-Likeness Prediction
Drug-like prediction of the compound was done using online SWISS-ADME [18]. Several different drug-like rules were evaluated for the compound like Lipinski rule of five, Ghose rule, Veber rule, Egan rule, and Muegge rule.

ADMET Property Prediction
ADMET stands for absorption, distribution, metabolism, excretion, and toxicity and forms pharmacokinetic profile of a given compound [19]. To predict these properties of the compound, an online admetSAR prediction tool was used [20].

Hirshfeld Surface Analysis
The conventional mapping of d norm and 2D fingerprint plots of newly synthesized three different heterocyclic compounds was drawn by Crystal Explorer 17.5 software with the aid of Crystallographic Information File (CIF). And, the cluster of molecules within the 3.8 Å radius were selected around a reference molecule to calculate pairwise interaction energies in a crystal, and energy frameworks were also carried out using CE-B3LYP/6-31G(d,p) method.
In the d norm surface, the strong red surface observed on nitrile and hydroxyl groups confirms O-H‧‧‧N type hydrogen bonding interactions in the crystal phase. The small red spots are also detected around the molecules; it indicates the C-H···N type of weak hydrogen bonding interactions. The strong white and blue surfaces distinguish the sum of van der Waals radii and short contact distances on the 3D d norm -based Hirshfeld surface maps. The shaped index shows the blue and red color which reflects the cycle stacking (C − H•••π and π•••π interactions). The two-dimensional fingerprint plots were created for each surface point and the data obtained from discrete intervals of d i and d e . In which, the different types of intermolecular interactions of the molecule in the crystal can be characterized by the shapes of fingerprint plots; the fingerprint plots of the titled compounds are shown in Fig (Fig. 2); the pairwise interaction energies were also estimated between molecules within a standard cluster of a radius of 3.8 Å at the B3LYP/6-31G(d,p) level of theory (CE-B3LYP model). Here, the total interaction energy between any nearestneighbor molecular pairs is given in terms of four components: electrostatic, polarization, dispersion, and exchange-repulsion which gives an insight into the underlying interaction energy and leads to the organization of crystal packing into supramolecular architectures in crystalline materials, which are summarized in Table 2. The representative energy framework diagram of pairwise interaction energies is shown as cylinders between molecular centroids; the thickness of each cylinder is proportional to the relative strength of interaction energies. The total energies of all interacting molecules with respect to the corresponding reference molecule along with the symmetry (x, y, z : − 39.2; − x, − y, − z: − 30.0; and − x, y + 1/2, − z + 1/2: − 41.8 kcal/mol) carry the highest energy than the other molecules and centroid-centroid distance. The electrostatic, polarization energies are lesser than the dispersion and exchange-repulsion energies, which reveals the topology of the molecules, and it concludes that these interacting energies are playing a major role in the assembly of the molecules in the crystal phase of molecules (Fig. 3).

Optimized Structure
The atomic charges, HOMO, LUMO, and electrostatic potential were also analyzed to understand the insight beauty of the molecule. And, the natural bonding orbitals (NBO) were also carried out to reveal the interactions that take a place between the filled and vacant orbitals. The graphical interpretation of all these studies was made by Gauss view and 3plot program. The comparison of optimized geometry of titled molecule reveals the similarity with crystal structure, shown in overlaid structure (Fig. 4). Table 2 Interaction energies (kJ/mol); R is the distance between molecular centroids (mean atomic posi- 1 -x, -y, -z 11.08 B3LYP

Global Reactivity Descriptors
To understand the relationship between structure, stability, and global chemical reactivity, the reactive descriptors of the titled compound can be derived from the conceptual density functional theory. The ionization potential and electron affinity are related to the energy of the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO). These orbital energies are highly helpful in a computational way using Koopman's theorem for closed-shell molecules to reproduce the reactivity of the molecule. The calculated value of ionization potential, electron affinity, and electronegativity of the titled compound are 7.39, 1.51, and 4.45 eV respectively. And, the calculated electrophilicity index (w) is 3.37 eV. The HOMO and LUMO are the main parts in the characterization of molecules in  (Table 3).

Atomic Charges and Electrostatic Potential
The knowledge of charge distribution of the molecule is very much essential in the field of quantum chemistry to the molecular system due to the effect of atomic charges, dipole moment, polarizability, etc. The creation of acceptor and donor pairs placed by the distribution of atomic charges over the molecule which is also connecting the charge transfer in and between the molecules. This charge transfer helps us to comprehend the chemical reactivity, molecular electrostatic potential, and electrostatic interactions. Atomic charges of the molecule are calculated to explore the electrostatic properties of the molecule which is calculated from the spherical charge approximation. The atomic charges of carbon and nitrogen atoms in the C≡N bond are approximately 0.307e and − 0.269e respectively. Among all the carbon atoms, methyl group carbon atoms carried a high negative charge (~ − 0.607e), and oxygen bonded keto carbon (C = O) carries high positive charges (0.843e). The more positive charges of the carbon atom are confirmed due to the substitution of the negative charge of the oxygen and nitrogen atoms in the molecule. Furthermore, the hydroxyl group  (O-H) oxygen atom possesses high negative charge (− 0.761e) than the other two oxygen atoms in the molecule, and the hydroxyl group hydrogen atom exhibits a higher charge than the other hydrogen atoms in the molecule. In this molecule, the hydroxyl group plays a major role in the inter and intra-molecular interactions. The molecular electrostatic potential is an important property to understand the molecular interactions and their charge distribution. Interestingly, it is very helpful to interpret and also predict the sites of electrophilic and nucleophilic attack. Therefore, the molecular electrostatic potential map plays a vital role in investigating the molecular structure along with the physiochemical property relationship. Figure 6 displays the electrostatic potential (ESP) map of the titled molecule; it also provides information about the role in drug-receptor interactions and the alignment of molecule in the receptor active site. A large negative ESP (red surface) corresponds to the attraction of proton by the concentrated electron density in the molecule, which is found around the C≡N bond region and C = O bond region in the molecule, and it acts as a nucleophilic attack. Whereas, positive ESP (blue surface) corresponds to the repulsion of the proton by the atomic nuclei and low electron density observed in these regions. Therefore, the molecular electrostatic potential map of the titled molecule confirms the negative surface around nitrogen and oxygen atoms while the positive surface around all carbon and hydrogen atoms. This may be used to predict the alignment of the molecule in the receptor active site.

Non-linear Optical Properties (NLO)
The organic molecules with second-order non-linear optical (NLO) properties exhibit a keen interest due to their significant role in the multidiscipline research like medicine, molecular switches, telecommunication, and data storage. There are many reports that confirmed certain electronic properties to have the desired nonlinear optical properties. In the computational way, the gaussian software gives different components of 3 × 3 × 3 matrix from the dipole moment (µ) of the molecule. Here, polarizability and hyperpolarizability of the titled molecule were calculated from the DFT/B3LYP 6-311G** method ( Table 4). The dipole moment, polarizability, small HOMO_LUMO energy gap, and planarity of the molecule are responsible for the non-linear optical properties. The dipole moment (µ), the (a) (b) Fig. 6 a ESP map and b atomic charges of the molecule polarizability (α total ), and first-order hyperpolarizability (β total ) of the molecule are 6.287 Debye, − 28.54 × 10 −24 esu, and 2.07 × 10 -30 esu respectively. Generally, the computed NLO properties of the organic materials are compared with urea which is also calculated at the same level of basis set; the value is 0.37 × 10 −30 . The hyperpolarizability (β total ) of the titled compound is much higher than the urea which clearly confirms the title molecule is a good material for the NLO applications.

Natural Bond Orbital Analysis (NBO)
The NBO analysis gives information about the inter-and intra-molecular bonding interactions to investigate hyperconjugation or delocalization interactions or charge transfers in the molecular system (Table 5). And, it also gives the details of the interactions between lone-pair-filled orbitals (bonding) and Rydgberg empty orbitals (anti-bonding).
The large E (2) value confirms the interaction between electron donors and electron acceptors and greater the extent of the conjugation of the whole system. The stabilizing donor-acceptor interactions are depending on the delocalization of electron density between occupied Lewis type and unoccupied non-Lewis type orbitals. The hyperconjugative interactions are the most responsible for the stability of the molecular system. The stabilization energies associated with electron delocalization between donor and acceptor are used to understand the inter-and intra-molecular interactions in the desired molecular system. From the NBO analysis, the titled compound forms three hydrogen bonds between the non-bonding orbital of the hydrogen-bonded acceptor (n A ) and the anti-bonding orbital of the H-donor bond (σ H−D* ). The stabilization energy of all three hydrogen bonds is greater than 0.5 kcal/mol.

Molecular Docking
The molecular docking technique is highly useful in predicting intermolecular binding conformation between molecules. The docking score of the compound is − 12.14 kcal/ mol and binds to the minor groove of the DNA. The binding mode of the compound 6.287 D Table 5 NBO analysis of quinoxaline derivative The type of n is as such that it positioned itself at the center of the minor groove; the ethyl formate moiety is more towards the leader while the central 2-hydroxy-2-methylcyclohexane-1,1-dicarbonitrile ring is placed straight over the nitrogen bases. The interactions of the compounds are dominated by van der Waals. The docked binding mode of the compound with DNA and its interactions are presented in Fig. 7.

Molecular Dynamics Simulations and Binding Free Energies
Molecular dynamics simulation of long 100 ns run was conducted to determine biding interactions and binding mode stability of the receptor DNA molecule in the presence of the drug. It has been unveiled that both DNA and the compound attained considerable stability along the simulation time as can be seen in Fig. 8. Conformation stability of the DNA and the compounds was evaluated via the root mean square deviations (RMSD) analysis.
The RMSD for the DNA and the compound fluctuates around 1 Å, which is an indication of a highly equilibrium system. The RMSD of the compound near 50 ns experienced a small deviation which corresponds to a bit conformational change in order to get more stability. Afterward, the RMSD plot for the compound was seen uniform and the binding conformation rigid. The DNA molecule in the presence of the compound enjoyed an enhanced level of stability. This data reflects on the high intermolecular strength and formation of a strong and stable complex. Additionally, to validate the predictions made by docking and simulation studies, MMGBSA binding free energies were estimated for the complex. As witnessed in the docking studies that the van der Waals energy contributes to the bulk of the interactions, the same was revealed in MMGBSA binding free energy analysis. The net binding van der Waals binding free energy of the complex is − 31.54 kcal/ mol while electrostatic energy is − 10.84 kcal/mol. The presence of electronegative atoms in the compound structure also allows the compound to produce electrostatic energy. Both the van der Waals and electrostatic energy is contributed high to the net binding energy. The net solvation energy of the system is 15.56 kcal/mol with favorable energy from nonpolar energy (− 9.48 kcal/mol) and non-favorable contribution from polar solvation energy (25.04 kcal/mol). The total binding energy of the system estimated is − 26.82 kcal/mol. The different binding free energy estimated for the complex is tabulated in Table 6.

Prediction of Drug Likeness
For a good drug-like molecule, it is necessary to be absorbed well in the host body in a given time. This supports good distribution of the compound and reaching the target site in maximum concentration for effective action [21]. Also, the drug needs to be  non-toxic as the failure of a drug molecule in clinical trials is highly costly. In silico prediction of such properties allow the selection of suitable drug molecule which does not accelerate the drug discovery process [22]. The compound was unveiled to fulfill all parameters of the Lipinski rule of five, Ghose rule, and Muegge rule. However, it fails by violating the TPSA value in both Veber and Egan rules. Each of the above rules uses its set of criteria for classifying whether a compound has a good chance to be a good drug candidate or not. According to the famous Lipinski rule of five, the compound TPSA value is within the required range and thus has better cell permeating ability and reaches the target site of action. The molecular weight of the compound is 438.48 g/mol. This weight affirms that the compound can get enhanced absorption rate [23]. The number of rotatable bonds, hydrogen bond acceptors, hydrogen bond donors, molar refractivity, and topological polar surface area (TPSA) is 5, 7, 1, 117.23, and 141.69 Å 2 , respectively. The LogP value of the compound is 2.83, which indicates good absorption and permeation of the compound. The compound has TPSA value within the limits and would be not a substrate for p-glycoprotein for efflux from the cell [24]. Details about fulfilling and violating a particular drug-like rules parameter are given in Table 7.

ADMET Property Analysis
The compound was unveiled to have better intestinal absorption with a probability score of 0.99. This demonstrates good intestinal absorption of the compound after oral administration. Similarly, the compound is predicted to have good blood-brain barrier (BBB) penetration. The compound is a non-inhibitor of P-glycoprotein (P-gp) and acted as a substrate of P-gp. From toxicity, the perspective of the compound is mutagenic, and hepatotoxic. Also, the compound showed oral toxicity. Furthermore, the compound was found to show zero alerts for PAINS thus can be ensured to show specific interactions with the DNA target only.

Conclusions
In summary, ethyl 3,3,5,5-tetracyano-2-hydroxy-2-methyl-4,6-diphenylcyclohexane-1-carboxylate was synthesized by three-component condensation of benzaldehyde with ethyl acetoacetate and malononitrile in the presence of trichloroacetic acid and its structure was proved by X-ray analysis. Then, several computational binding analyses have been carried out. Firstly, the mapped Hirshfeld surface analysis confirms that the O-H…N type of hydrogen bonding interactions is in the crystal phase. The methyl group carbon atoms and the hydroxyl group (O-H) oxygen atom exhibit high charges. The nitrile and keto groups were found to be electronegative regions. Furthermore, the NLO and NBO analyses of the molecule reveal more detailed molecular properties of the titled molecule. Computer-aided drug analysis demonstrated the good binding of the compound to the DNA molecule. Dynamically, the compound binding conformation with the DNA remained stable and interacted with good van der Waals energy. The compound was also reported to possess good drug-like properties and had a good pharmacokinetic profile; however, the compound contains toxic chemical moieties that need to optimize in the future.
Author Contribution Malahat Kurbanova, Arzu Sadigova, and Abel Magerramov participated in the emergence of the synthesis. Rizvan Askerov contributed to X-ray analysis and Youness El Bakri, Kandasamy Saravanan, and Sajjad Ahmad took part in theoretical calculations.

Availability of Data and Materials
All data and materials that were used in the article are available.

Declarations
Ethical Approval All ethical norms were maintained by the authors during the article preparation process.

Consent to Participate
All authors agree to participate in the process.

Consent for Publication
All authors agree to publish the article.

Competing Interests
The authors declare no competing interests.