DFT calculation, molecular docking, and molecular dynamics simulation study on substituted phenylacetamide and benzohydrazide derivatives

The atomic and molecular properties of the title compounds were calculated by Jaguar using a basis set B3LYP/6-31G**++ with hybrid DFT in the gas phase, to determine the chemical reactivity. Analysis of quantum chemical features such as HOMO and LUMO explained that the electronic charge transfer occurred within the system through conjugated paths of the selected compounds. The nucleophilic and electrophilic reactive sites are recognized from the molecular electrostatic potential plot. Electrophilic and nucleophilic attack-prone molecular sites were predicted by mapping ALIE value to the molecular surface. The bond dissociation energy of the high active compound 15 (2-chloro-N-(2-(2-(2-(2-chlorobenzoyl)hydrazineyl)-2-oxoethoxy)phenyl)acetamide) was calculated to assess the probability of compound autoxidation or degradation. Further, molecular docking, binding free energy calculations, and ADMET profile of the degradation products (DPs) of compound 15 was carried out to determine the binding affinity and toxicity profile of the formed DPs compared with the parent compound. A 150-ns molecular dynamics (MD) simulation was performed to evaluate the binding stability of the compound 15/4URL complex using Desmond. Binding free energy and binding affinity of the complex were computed for 100 trajectory frames using the MM-GBSA approach.


Introduction
Antibacterial resistance (ABR) is a rising global burden, threatening our ability to treat common infectious diseases [1,2]. As per the statistics by the center for disease control and prevention report 2019, 15 more antibiotic-resistant strains were discovered and classified as "urgent," "serious," and "concerning" threats. This demands the discovery of small-molecule inhibitors to overwhelm ABR. Consequently, various strategies have emerged to develop new molecules that could either act as drugs or by adjoining the different bioactive moieties that are known to possess biological activity against nosocomial pathogens [3][4][5][6][7].
In the recent past, phenylacetamides [8,9] and benzohydrazides [10][11][12] have garnered attention among researchers, especially in the discovery and synthesis of novel compounds with notable biological activities. These moieties contain active amide and hydrazide fragments that act as links between any aryl or heteroaryl groups. Due to their intensive pharmacological activity, these are widely used in organic synthesis. Phenylacetamides [13,14] and benzohydrazides [15] find application in designing and developing novel drugs that are known to possess antitubercular, antimicrobial, antimalarial [12], and antiinflammatory [16] properties. Phenylacetamides and benzohydrazides with electron-withdrawing and donating groups provide a pathway for redistribution of electron density under the influence of the external field. Such compounds are of utmost importance in material chemistry. Their strong chemical stability serves as an attractive future material in different areas such as optoelectronics technologies, optical switching device, and chemo sensor development.
Pharmaceutically, organic molecules are the candidates that are synthesized to be stable in the aquatic environment. Regrettably, a typical environment is not appropriate for the degradation of the organic compounds in water, which led to their accretion in all types of water worldwide [17][18][19]. Hence, the stability and reactivity of the molecules command the fate of the pharmaceutical ingredient in the environment [20].
However, forced degradation methods are performed to enhance the degradation techniques of active ingredients in the formulation [21]. Moreover, these studies are of great importance from the aspect of desirable resources; thus, the computational tools are employed for optimization and rationalization of such experiments [20,22,23]. Density functional theory (DFT), Ab initio, Semi-empirical Self-consistent Field (SCF), and Molecular mechanics approaches serve as the crucial tools of computational chemistry. In addition, DFT was used to calculate model structures such as fluorinated ethylene propylene (Teflon FEP), polymethyl methacrylate, and Kapton polyimiide (Fig. 1).
Further, specific correlations among different components are readily computed by in silico approaches from one side and the physicochemical properties from the other aspect. One example is bond dissociation energies (BDE) of DFT calculations which can be employed for the molecule site prediction where autoxidation might occur. The degradation products (DPs) formed during storage are the latent source for toxicity, which plays a crucial role in drug discovery and development. Thus, the characterization and toxicity evaluation of predicted DPs may provide support to lessen the adverse effects and avoid the recall of drugs from the market [24].
Moreover, quantum chemical calculations are reasonably necessary to understand electronic polarization, responsible for predicting novel properties and non-linearity. The non-linearity of the compound can be synthetically modeled by extending the π-conjugated system or by optimizing donor-acceptor strength. The effect of water and hydrolysis of tested molecules can be evaluated by molecular dynamic (MD) simulations and radial distribution function calculations.
Thus, considering the importance of phenylacetamides and benzohydrazides from the pharmaceutical aspect, we employed DFT calculations and MD simulations to anticipate some significant reactive properties. In our earlier studies, we reported the synthesis and antibacterial activities of phenylacetamides and benzohydrazides [25,26] (Supplementary Table S1). Here we report the DFT calculations and MD simulation studies on the above-mentioned molecules in continuation of the above work.

DFT calculations and electrostatic potential (ESP) calculations
All atomic and molecular properties of the compounds were calculated by Jaguar (v8.8, Schrödinger 2019-2) using a basis set 6-31G** ++ with hybrid DFT of Hamiltonian-nonrelativistic correlation function (B3LYP) in the gas phase (Table 1). DFT energetics demonstrates 3D electronic states of molecules to evaluate the transfer of lone pairs, bonds, and reactivity in a specific environment [27]. Further, selected ligands were examined for their quantum chemical features, such as localization energies of lowest unoccupied molecular orbital (LUMO) and highest occupied molecular orbital (HOMO) and molecular electrostatic potential (MESP) using Jaguar [28]. These localization energies are named frontier molecular orbitals (FMOs) that display a crucial role in maintaining chemical stability and are also used as an efficient tool for investigating donor-acceptor interactions [29]. HOMO characterizes the ligand's ability to donate an electron, whereas LUMO refers to the ability of the ligand to accept an electron. Low values of LUMO indicate the higher propensity of the molecule to accept electrons, whereas higher values of HOMO dictate the increased propensity to donate electrons to unoccupied molecular orbitals. The values of HOMO, LUMO, and their energy gap reproduce the level of conductivity, reactivity, and kinetic stability of the compound. In soft molecules, electron transfer from HOMO to LUMO with low kinetic stability and high chemical reactivity is more effortless if the value of the energy gap is small, which influences the biological activity. The higher energy gap between HOMO and LUMO indicates the kinetic instability of molecule, whereas narrow HOMO and LUMO energy gap enables intra-molecular charge transfer, making the compound's non-linear optical (NLO) properties active [30,31]. Structural geometry of ligand was optimized at B3LYP/6-31G** ++ level, and energy calculation was performed with Poisson-Boltzmann Finite (PBF) solvation state [32]. The DFT energy of a system is calculated as where E NN is the nuclear-nuclear repulsion, E COUL is the electron-electron Coulomb repulsion, E V is the nuclear electron attraction, E ET is the kinetic energy, E CORR is the correlated movement of electrons with different spin, and E EXCH is the electron-electron exchange energy.
MESP demonstrates the electron density distribution around the molecule and also predicts the molecule chemical reactivity. MESP is a tool for predicting chemical reactivity where negative potential regions (protonation sites) involve in the nucleophilic attack, whereas positive potential region is susceptible for electrophilic attack. DFT energy and MESP knowledge are mainly employed to predict possible functional spots in ligands for potential interactions with the residues of the target. Moreover, potential molecules were evaluated in terms of the HOMO and LUMO energy gap [33].

Bond dissociation energy
Further, bond dissociation energy (BDE) of the high active molecule was performed using Materials Science suite (Schrödinger 2019-2). Geometry optimization was performed without considering any constraints using hybrid density functional B3LYP with 6-31G** ++ basis set.

Molecular docking and binding free energy calculations
Docking studies were performed in extra-precision (XP) mode [34] for compound 15 and its probable degradation products (DPs)/BDE fragments using Glide module (Schrödinger 2019-2) without applying any constraints. The crystal structure of Staphylococcus aureus ParE enzyme complexed with keibdelomycin (PDB ID: 4URL, resolution 1.9 Å) was retrieved from the RCSB-PDB and then prepared using the Protein Preparation Wizard module. Ionization states of amino acid residues at pH 7.0 were generated, and hydrogen atoms were added to the protein. Missing side chains or loops were filled using the module Prime [35]. Energy minimization and optimization were performed using the OPLS3e force field [36] with root-mean-square deviation (RMSD) at 0.30 Å. In this study, a glide grid was created at the centroid of the co-crystal ligand keeping the partial charge cut-off at 0.25 and van der Waals scaling of 1.0 for the receptor. Compound 15 and its DPs were optimized using the LigPrep module to generate low-energy conformers. Binding free energy calculations of selected ligands were calculated using the Prime, MM-GBSA approach. The energy calculations of protein-ligand complexes were performed using the continuum solvent (VSGB 2.0) model [37] and OPLS3e force field.

Prediction of pharmacological and toxicological properties
In recent years, the significance of favorable pharmacokinetic features for discovering successful drug candidates has been widely acknowledged. Thus, ADMET evaluations are integrated earlier into the field of drug design and discovery. QikProp [38] proved to be a novel tool in computing and optimizing the pharmacokinetic profile or ADMET relevant descriptors for pharmaceutically acceptable compounds. Further, this algorithm is also used to correlate physicochemical and pharmacokinetic properties with the 3D molecular structure of the test ligands. Neutralization of compounds is an essential step that is performed before being used by QikProp, as this tool is incapable of neutralizing a compound. A total of 44 predicted properties such as physicochemical properties, principal descriptors, LogP, QP%, and log HERG are generated in normal mode. It also determines the acceptability of the ligands, based on Lipinski's rule of five [39] and three, which are crucial for rational drug design.

Molecular dynamics simulation and analysis
MD simulations were performed of the docked complex 15/4URL using the OPLS3e force field [35]. The system was solvated in an orthorhombic boundary box with explicit TIP4P water within the Desmond MD system [40,41]. Counter ions were added to achieve the neutralization of the system. The system comprises 50,390 atoms and 14,800 water molecules. System minimization was carried out with OPLS3e force field using 200 steepest descent algorithms, followed by 1000 conjugate gradient steps until a gradient threshold of 25 kcal/mol reached. The Smooth Particle Mesh Ewald method was employed for short-range van der Waals and Coulomb interactions at a cut-off radius of 9 Å, while tolerance of 1e − 09 was applied for the long-range electrostatic interactions [42]. A 150-ns MD simulation was performed under isobaric and isothermal ensemble (NPT) with a pressure of 1 bar and temperature of 300 K using Martyna-Tobias Klein and Noose-Hoover thermostat methods, respectively [43,44]. RESPA integration algorithm [45] was employed for the bonded, near, and far non-bonded interactions with multiple time steps of 2, 2, and 6 fs, respectively. Maestro graphical user interface was used to visualize the trajectories and 3D structures. Further, binding free energy calculation of 15/4URL complex was carried out using the MM-GBSA approach.

Results and discussion
The atomic and molecular properties of the title compounds were calculated using Jaguar (Schrödinger 2019-2) using a basis set 6-31G** ++ with hybrid DFT with Hamiltoniannonrelativistic correlation function (B3LYP) in the gas phase, and the results are summarized in Table 1. In this article we discussed DFT calculations, bond dissociation energies mainly focusing on compound 15 which considered to be highly active among the tested compounds.

Frontier molecular orbital analysis
The frontier molecular orbital analysis demonstrates the intramolecular charge transfer between the acceptor and donor functional groups through the conjugated paths [46]. In the current study the HOMO and LUMO energies are in the range of 9.57319 to 4.58334 eV and 2.19099 to 0.12709 eV, respectively. The HOMO of the high active compound 15 (MIC, 0.62 µg/ml) is localized at − 9.060094 eV, while the LUMO is localized at − 5.585780 eV, and there is a charge transfer in the molecular system through the conjugated paths. The energy gap (ΔE) of − 3.474314 eV between HOMO and LUMO specifies the molecular chemical stability (Fig. 2). Moreover, the ΔE values showed good correlation (R 2 = 0.812) with the MIC of compounds 1-29 (0.68-8.99 µg/ml) against Staphylococcus aureus (NCIM 5022).

Molecular properties
Electronegativity is a measure of the propensity of an atom to attract electrons in a chemical bond or negative chemical potential of an atom. In contrast, chemical hardness is a measure of resistance to charge transfer. The electrophilicity index of a system demonstrates the energy stabilization and electron affinity and measures the binding energy decay resulting from electron flow between a donor and an acceptor. The molecular electronegativity and hardness of the selected compounds were found to be in the range of − 5.601 to − 1.586 a.u. and 4.3256 to 2.699 a.u., respectively (Table 1). However, the high active compound 15 exhibited maximum value of electronegativity (− 4.802 a.u.) and hardness (4.234 eV).
Moreover, the reactivity of the carbonyl group of a molecule can be explained by the total nucleophilic super delocalizability (TNS) of the carbonyl group, whereas reactivity of olefinic double bond is measured by electrophilic total super delocalizability (TES) [47]. From Table 1, it is evident that the TNS (TNS, − 131.326) is more significant than TES (− 7.454) of olefinic double bond. This implies that carbonyl group of compound 15 behaves as a nucleophile rather than an electrophile [47].
Molecular electrostatic potential (MESP) is an essential property for predicting and analyzing the molecular reactivity behavior [48] of selected ligands. MESP was computed for compound 15 (Table 1) to indicate the relative reactivity sites for the nucleophilic or electrophilic attack, hydrogen bonding interactions, crystal behavior, and macroscopic properties. The positive (light blue) and negative (red) regions in the MESP plot (Fig. 3a) indicates nucleophilic and electrophilic reactivity, respectively. The high active compound 15 exhibited several possible sites of both electrophilic and nucleophilic attack. However, negative ESP points were mainly localized over the carbonyl oxygen and electron cloud of the phenyl ring and the value was found to be 0.07 a.u. on the molecular surface. The deepest positive ESP is localized on the nitrogen and hydrogen atoms with + 0.52 a.u. All these sites indicated the possible strong intermolecular and intramolecular interactions in the crystal lattice. Further, identifying molecular sites prone to electrophilic attacks is of utmost importance for the overall understanding of the molecule's reactivity. ALIE is the energy needed to remove an electron from a specific point or in the space of the system. Average local ionization energy (ALIE) surface was generated [49] to determine the molecular sites where electrons are loosely bound and tightly bound. Analysis of ALIE of the high active compound 15 explained that carbonyl oxygen atoms are the locations of the loosely bound electrons and hence favorable for the electrophilic reaction or radical attack (Fig. 3b).

Bond dissociation energy (BDE) of compound 15
For the higher BDE of compound more energy is required for the bond cleavage and vice versa. A total of 16 DPs were formed for compound 15 by 8 dissociation patterns. The results show that compound 15 exhibited low BDE for DP1a and DP1b (Fig. 4), and the value was found to be 46.161 kcal/mol. In contrast, a high BDE of 85.406 kcal/mol was observed for the DP8a and DP8b.   Figure S1). The compound 15 exhibited three hydrogen bonding interactions with Glu53, Gly80, and Thr168. Additionally, this compound was stabilized by π-π stacking interaction of phenyl ring linked to the acetamide fragment with the phenyl ring of Phe106. Further, a π-cation interaction was observed between the same phenyl ring of compound with the protonated NH 2 of Arg79 (Fig. 5). The glide scores for DPs were observed to be in the range of − 5.011 to 4.001 kcal/mol. Prime-based MM-GBSA result indicated high binding affinity within the key binding pocket residues of ParE enzyme. The calculated ΔG bind values were found to be in the range of − 63.51 to − 12.44 kcal/mol. Further, hydrophobic energy (ΔG lipo : − 19.88 to − 7.65 kcal/mol) term was found to be favorable for the binding of the ligands with ParE protein. Coulomb (ΔG cou : − 17.03 to 34.44 kcal/ mol) and covalent energy (ΔG cov : − 0.60 to 19.86 kcal/mol) terms were observed to be moderately unfavorable for the binding.

Pharmacological and toxicological parameters
The ADMET result obtained through QikProp is represented in supplementary Table S3. ADMET properties of compound 15 and its DPs were found to be in the recommended range. The analysis of the results demonstrated the similarity in the ADMET profile of compound 15 and its DPs. The results exhibited that compound 15 and its DPs possess druggable properties and obeyed the Lipinski rule of five with 0 to 1 violation. The predicted CNS activity of compound 15 and its DPs was found to be in the range of (+ 2 as active and − 2 as inactive), exmplifying a lack of CNS effect. Partition coefficient (QPlog O/W ) of compound 15 and its DPS (0.2020 to 3.498) is well within the recommended range. Besides, the predicted apparent Caco-2 cell permeability of compound 15 and its DPs was found to be in the range of 162.774 to 1449.673 nm/s, indicating the high active permeability of these compounds. Solvent accessible surface area (SASA) of compound 15 and its DPs were also observed to be in the recommended range (357.577 to 699.599 Å 2 ). Compound 15 possesses greater SASA (699.599 Å 2 ), indicating the burial of this compound within the binding pocket compared to its DPs. Compound 15 and its DPs did not exhibit blockadge of HERG K + (− 2.964 to − 3.409) indicating the safety of this compound. Moreover, the polar surface area (PSA) of compound 15 and its DPs were found to be in the range of 56.161 to 120.338 Å 2 which is within the recommended range (7-200 Å 2 ), which shows better van der Waal surface area around these molecules.
Further, compound 15 and its DPs possess greater topological surface area (TSA) within the recommended range (< 120 Å 2 ), which indicates the total sum of overall polar atoms or molecules (primarily oxygen and nitrogen, including their attached hydrogen atoms). The compound 15 and its DPs exhibited non-permeability through the blood-brain barrier and are not substrates for P-gp. Further, these are not considered as cytochrome P450. LogKp is a skin permeation indicator showing the compound absorption through the skin. The compound 15 and its DPs possess good permeation through the skin (− 6.58 to − 7.88 cm/s). However, all DPs possess a bioavailability score (0.55) similar to that of the parent compound. Additionally, compounds exhibited zero Pains alert (PA) (PAN assay interference compounds), indicating the specific interaction with desired target and no interaction with other biological targets.

Molecular dynamics simulation and analysis
The molecular flexibility, structural behavior, and stability of ParE docked with high active compound 15 was assessed by 150 ns of MD simulation using Desmond. The root-meansquare deviation (RMSD) values of protein Cα backbone and heavy atoms increased sharply up to 36 ns and then stabilized in the narrow range of 2.53-3.90 Å, 2.54-3.90 Å, and 2.98-4.09 Å, respectively, during rest of the simulation time (Fig. 6a). This indicated less fluctuations in protein structure during MD simulation.
The ligand root-mean-square fluctuation (RMSF) values of Cα atoms, backbone, and heavy atoms of catalytic pocket binding residues were found to be in the range of 0.64-3.01 Å, 0.65-2.90 Å, and 0.76-2.90 Å (Fig. 6b). The maximum RMSF was observed for Phe106 and Lys112 from N-terminal domain present in the loop ranging from Val99 to Gly121, connecting two α-helices Val122 to Ala127 and Val93 to Thr98. The residues contacting to ligand were Fig. 6 a Root-mean-square deviation (Å) of Cα atoms and backbone. b Root-mean-square fluctuation profile. c Interaction profile of complex 15/4URL during 150-ns simulation observed between Ile46 and Val170 due to pronounce effect of catalytic activity in this region. Whereas no contacts were found in the region between Thr171 and Arg399.
Further, the MD trajectory of 15/4URL protein complex revealed hydrogen bonding interactions mainly with the polar and charged residues of N-terminal domain dominated by hydrophilic residues. Among the four hydrogen bonds observed in XP-docking, only one (Glu53) was preserved during MD simulations. This may be due to the evolution of inhibitor during MD simulation. Interactions were observed with residues stretching between Ile46 to Val170 (Figs. 6c and 7). However, no π-π stacking, π-cation, and ionic interactions were observed during 150-ns MD simulation. A total of 8 hydrogen bonding interactions were observed, out of which Gly104 formed two strong hydrogen bonds with the nitrogen of phenylacetamide (>C=O⋯HN, 77% of MD trajectory) and the carbonyl oxygen of hydrazide linker (HN⋯O=C< , 76% of MD trajectory). Precisely, carbonyl moiety of phenylacetamide formed a strong hydrogen bonding interaction with Thr92 via a water bridge (>C= O⋯H⋯O(H)⋯HN, 66% of MD trajectory). The two NH of hydrazide linker established low-frequency hydrogen bonding interaction with Asp52 (18% of MD trajectory) and a moderate-frequency hydrogen bonding with Asn49 (38% of MD trajectory) via water bridge. In addition, the carbonyl group of hydrazide linker accepted a water-mediated hydrogen bonding interaction with Asp76 (>C=O⋯H⋯O(H)⋯HN, 22% of MD trajectory). Additionally, ligand also exhibited low-frequency hydrogen bonding interactions with Glu53, Asp76, Arg79, Met81, Pro82, Ile96, and Phe106 (17-19% of the MD trajectory). However, π-π stacking and π-cation interactions were not observed for this compound. It is evident from Supplementary Figure S2 that the phenyl ring of the phenoxy moiety is buried between hydrophobic residues, whereas hydrazide linker and polar fragment (−NH-C=O-CH 2 -Cl) present at position two of the phenyl ring of the phenoxy moiety is completely embedded within the hydrophilic residues of the catalytic pocket. It is also evident from the above result that Gly104, Thr92, Asn49, and Asp52 played a crucial role in stabilizing compound 15 within the binding pocket of 4URL. Further, XP-docking pose of compound 15 was overlayed with the MD simulation pose of this compound with an RMSD of 1.23 Å (Supplementary Figure S3). This validates our docking protocol and show the stability of this compound within the catalytic pocket.
The radius of gyration is indicative of the compaction level, or extendedness of the ligand was found in the range  Figure S5).
We computed binding free energy from the generated 100 trajectory frames of 150-ns MD simulation using the module thermal_mmgbsa.py Script. The average ΔG bind was observed to be − 60.35 kcal/mol, which correlates well with the computed MM-GBSA (− 63.51 kcal/mol) (Supplementary Figure S6).

Conclusion
HOMO and LUMO values of compound 15 explained that the electronic charge transfer occurred within the system through conjugated paths. The ΔE HOMO-LUMO of compound 15 was found to be − 3.474314 kcal/mol, indicating the kinetic stability of the molecule. Further, MESP studies revealed that the carbonyl oxygen of hydrazide moiety and oxygen atom of phenol ring are the most negative electrostatic potential regions, while nitrogen of the − NH − fragment possess positive electrostatic potential. BDE values represent possible degradation products whose values are greater than 45 kcal/mol. Further, the toxicity profile of degradation products was compared with the compound 15. Based on the results obtained from docking and free energy studies, it is evident that DPs have less affinity for the binding residues of catalytic pocket compared to the compound 15. A 150 ns of MD simulation was performed to validate the binding stability of the 15/4URL complex. From the results, it is evident that the stability of the ligand was mainly dependent on the interactions with polar and charged residues Gly104, Thr92, Asn49, and Asp52 from the N-terminal domain.