Validation of crystal structure of 2‐acetamidophenyl acetate: an experimental and theoretical study

Abstract In this present study, we have determined the crystal structure of 2-acetamidophenyl acetate (2-AAPA) commonly used as influenza neuraminidase inhibitor, to analyze the polymorphism. Molecular docking and molecular dynamics have been performed for the 2-AAPA-neuraminidase complex as the ester-derived benzoic group shows several biological properties. The X-ray diffraction studies confirmed that the 2-AAPA crystals are stabilized by N–H···O type of intermolecular interactions. Possible conformers of 2-AAPA crystal structures were computationally predicted by ab initio methods and the stable crystal structure was identified. Hirshfeld surface analysis of both experimental and predicted crystal structure exhibits the intermolecular interactions associated with 2D fingerprint plots. The lowest docking score and intermolecular interactions of 2-AAPA molecule against influenza neuraminidase confirm the binding affinity of the 2-AAPA crystals. The quantum theory of atoms in molecules analysis of these intermolecular interactions was implemented to understand the charge density redistribution of the molecule in the active site of influenza neuraminidase to validate the strength of the interactions. Communicated by Ramaswamy H. Sarma


Introduction
Influenza is an acute respiratory disease that affects between 20 and 40 million people per year in the United States, mostly the immunocompromised population. The 1918 (H1N1 "Spanish"), 1957 (H2N2 "Asian"), and 1968 (H3N2 "Hong Kong") influenza pandemics were all caused by the type A influenza virus (Hsieh et al., 2006). The influenza virus membrane has two glycoproteins; hemagglutinin (HA) and neuraminidase (NA), both recognize the sialic acid (AJ et al., 2001). Commencement of viral infection is due to the binding of different HAs with sialic acids on the cell surface. After the host cell got infected, NA, breaks the binding with sialic acid, and spreads the infection to other cells. Thus, the NA intercedes an essential advance in the existence pattern of influenza A virus, the focus of this paper is antiviral drug complexes with NA enzymes.
The NA enzyme is ordered into eleven subtypes (N1 to N11). They are gathered under two groups; group 1 N1, N4, N5 and N8 and group 2 contains N2, N3, N6, N7 and N9 structure (SJ & JJ, 2010). Both the groups were distinguished from the conformational differences in the active site loop segments, spotted by residues 148-151. In Group1, NAs exhibits the open conformation state and have space for ligand inhibition, whereas in the group 2, NAs pose closed conformation at the active sites.
M arquez-Dom ınguez et al. (2020) reviewed the active site of N1 complexed with the antiviral drugs zanamivir. The active sites of N1 were grouped into five subsites (S1-S5). The three essential amino acids (R118, R292, and R371), of S1 and amino acids; E119, D151, and E227 of S2 were considered as crucial segments engaged in hydrogen bond network in the active site. The interactions in other subsites are hydrophobic in nature.
The emergence of drug-resistant viruses requires the development of new antiviral chemicals. Benzoic acid derivatives have been reported to possess anti-influenza virus activities. In this view, the N1 enzyme interactions with 2 acetamidophenyl acetate [2-AAPA] (Scheme 1) was reported first time in this paper and compared the same with two reported drugs Zanamivir (Xu et al., 2008) and oseltamivir (Shi et al., 2006). However, to understand its behavior in the active site, detailed studies on the molecule by quantum chemical models and molecular dynamic simulation are essential.
2-acetamidophenyl acetate (2-AAPA) crystals were grown in ethanol by a slow evaporation method, among the harvested crystals, a good quality single crystal (ca. 0.780 Â 0.330 Â 0.090 mm) was chosen for a high-resolution X-ray diffraction experiment. Data collection of the title compound was performed using STOE IPDS 2 (Tanak et al., 2009) diffractometer provided with graphite monochromatic MoKa (k ¼ 0.71073 Å) radiation at 293 K. The intensities were corrected for Lorentz and polarization effects, and an analytical absorption correction was applied, using SCALE3 ABSPACK in CrysAlisPro (H€ ubschle et al., 2011). The structure was solved using SHELXS-97 (Sheldrick, 1997) and the data were refined by full-matrix least-squares methods on F 2 using the SHELXL-97 (Sheldrick, 1997) program package. The refinement details are given in Table 1.

Molecular docking and molecular dynamics simulation
To perform molecular docking, the 2-AAPA molecule was minimized with OPLS_2005 (Optimized Potentials for Liquid Simulations) (Petrovic et al., 2010) force field using the Ligprep application incorporate in the Schrodinger package (Wizard, 2014). The three-dimensional coordinates of the crystal structure of neuraminidase (PDB:3B7E) (Xu et al., 2008) were retrieved from a protein data bank. Further, to relieve steric clashes, restrained minimization was performed using a protein preparation (Wizard, 2014). To refine the conformational flexibility of protein while ligand binding, induced fit docking (IFD) was performed, which is integrated into Schr€ odinger software 8.0, LLC (Schr€ odinger, 2005). The best protein-ligand complex was identified based on the scoring and intermolecular interactions. The intermolecular interactions between the protein and ligand were analyzed from PyMol (Seeliger & de Groot, 2010). Using YASARA structure and dynamics, version 20.12.24.W.64 (Krieger et al., 2004(Krieger et al., , 2006, with the AMBER14 forcefield (Maier et al., 2015) molecular dynamics simulation was performed for the NA-2-AAPA complex. At a temperature of 298 K, the complex was solubilized in a cubic simulation cell that was 20 times larger on each side than the NA-2-AAPA complex. The water solvent density was 0.997 g/ ml. The boundary of the simulation cell is periodic. Based on the predicted pKa of the residues side chain from the Ewald summation, the charged residues were assigned (Krieger et al., 2006). About the physiological pH 7.4, hydrogen atoms were incorporated into the protein structure according to the computed pKa. To ensure the electroneutrality of the system, 41 Na þ ions and 37 Clions were added to the final concentration. Using the steepest descent minimization method, followed by simulated annealing, the structure was then minimized. For Lennard-Jones forces (Errington et al., 2003), a cutoff value of 8 Å was used while the Coulomb forces associated with the concerned structures were calculated using the Particle Mesh Ewald algorithm. The system was initially minimized using simulated annealing. After annealing, MD simulations were run at 298 K for 135 ns. For each 100 ps, the simulation snapshot for the trajectory was obtained and analysed using the YASARA structure package. Furthermore, free binding energies of NA-2-AAPA complex have been calculated at regular intervals of 10 ns with molecular mechanics-generalized born surface area (MM-GBSA) protocol using the prime module of Schrodinger software (Bhachoo & Beuming, 2017;Hayes, 2012).

X-ray crystal structure
The molecular structure of the title compound is shown in Figure 1. From the C-N bond lengths, C3-N1 [1.406(7) Å] and C4-N1 [1.358(7) Å], the electron delocalization is evident and the same was further proved from the C3-N1-C4 bond angle [126.6(5)˚] which is greater than the usual tetrahedral angle for nitrogen atom. The partial double-bond character of the C4-O3 [1.204(8)] and C1-O2 [1.192(7)] bonds are close in duration to those observed in previous reports [1.215(2)]. (Kansı z & Dege, 2018). The Phenyl ring lies almost perpendicular [86.9˚] with methyl acetate and 20.6˚with methyl acetamide groups. The root means square deviations of C-C-C bond angle in the phenyl ring from an ideal value (120 ) is calculated from r CÀCÀC ¼ ͱ P 6 i¼1 ðhiÀ120 Þ 2 6 , and the value is 1.44T he molecular conformation in the crystal is maintained by a strong intramolecular N-HÁÁÁO hydrogen bond. The intra and intermolecular hydrogen bonds prevailed in the system ( Figure 2) which stabilizes the crystal phase are tabulated in Table 2. The oxygen atom (O2) of the methyl acetate group forms strong hydrogen bond interactions with nitrogen atom [N1] of methyl acetamide and carbon atom [C7] of methyl acetate groups and forms R 1 2 (6), and R 1 2 (9) ring motifs. There is no pÁÁÁp and C-HÁÁÁp interactions between the molecules, the crystal structure was further stabilized by weak C-HÁÁÁO hydrogen bonds.

Conformers of 2-AAPA
The hypothetical stable conformers of 2-AAPA molecules generated from the DMACRYS optimization are tabulated in Table 3 The energy landscape scatters plot ( Figure S1), along with the comparison studies. indicated that theoretical structure FA, generated at the À96 KJ/mol, shows a strong agreement with the experimental crystal form of 2-AAPA. The similarity in the density (1.284/1.241) g/cm 3 ] and cell volume (978/1033) g/cm 3 further supports that structure FA is equivalent to the experimental 2-AAPA; authenticating the ab initio CSP methodology. The conformity of the theoretical structure FA with the experimental one was validated from the investigation of the possible hydrogen bonds associated with the crystal phase. Table 2 shows the possible hydrogen bonds associated with theoretical 2-AAPA crystal and those obtained experimentally. The predicted hydrogen bonds agree well with those of the actual 2-AAPA crystal, which is responsible for the stability. The variation of $0.5 Å between the donor-acceptor bond distances of the theoretically predicted molecule and the experimental values indicates the identical nature of the predicted structure with the experimental 2-AAPA. The studies indicated the hydrogen bond between nitrogen and oxygen atom in 2-AAPA (N-HÁÁÁO) as the crucial linkage that decides the stability of the crystal phase. The similarities in the lattice parameters for experimental and theoretical crystal forms of 2 AAPA show the authenticity of the methodology, hence validating the studies.

Hirshfeld surface analysis
The role of fundamental interatomic interactions in the stability of both experimental (real) and predicted (FA) crystal phases of 2-AAPA was confirmed from the 2D fingerprint plots generated using Crystal Explorer 17.5 software (Spackman & McKinnon, 2002). The term 'd norm ' leads to different factors: d e , which refers to the distance between any surface point to the interior of atoms, and d i , which refers to the distance between any surface point and exterior point of atoms, as well as the atom's van der Waals (vdW) radii. Interactions are outlined using a red-white-blue color code (shorter than vdWequivalent to vdWlonger than vdW isolation, respectively (Madura et al., 2010;Seth et al., 2011). The Hirshfeld surfaces mapped over d norm (range of À0.4077-1.4296 Å). The surfaces are shown as clear to permit the visual image of the molecule. The maximum contribution over the total Hirshfeld surface was encountered from HÁÁÁH contributes about 47% of the interactions for both experimental and predicted crystal structures. Next to this, the Hirshfeld surface (HS) was the pointed nature of the fingerprint plots indicated the dominant nature of strong OÁÁÁH/ ÁÁÁO interactions (Figure 3a,b). The corresponding contribution in percentage for experimental crystal structure is 31.3% and this value is slightly decreased to 26.8% for the  predicted crystal structure. Figure 3c displays the fingerprint plot difference between the experimental and predicted molecule. Blue regions are more intense in the experimental structure and red region highlights for the predicted one. The less dominancy of OÁÁÁH/ÁÁÁO contact in HS of the predicted structure was revealed near the spike at de ¼ 1. Figure 3c. The results from the contact enrichment ratio (Table 4) established a similar pattern in the interaction energies. Also, it is found that there is 96% correlation similarity between the intermolecular interaction for generated conformer and empirical structures of 2 AAPA ( Figure S2), which reveals the structures are equivalent. For both the experimental and projected crystal structures of the 2AAPA molecule, Figure 4a,b shows a scatter plot graph of the inner and outer electrostatic potential V on the Hirshfeld surface. The points on the surface are distinguished by interaction type. In Figure 4a, the region of maximum inner potential represents the polar hydrogen atoms, H O , and the oxygen atoms were at the highest exterior potential. A similar trend was found for the predicted crystal structure as shown in Figure 4b. The experimental and predicted crystal shows an electrostatic complementarity as a correlation C vv of À0.48% and À0.18%, respectively, were seen. The contacts of polar hydrogen atoms H c result in the most electropositive regions and represent strong hydrogen bonds OÁÁÁHn and NÁÁÁHn. The C-HÁÁÁO and C-HÁÁÁN contacts are weaker hydrogen bonds involving the most electronegative regions and mildly electropositive regions.

Interaction energies
From the pixel calculation, it is analyzed that the crystal packing of both experimental and predicted structure is    highly stabilized by dispersion energies, À102 and À94.1 kJ/ mol, respectively. The intermolecular energies between the molecular pairs are listed in Table 5. Among the energies of intermolecular interactions, the crystal stability was highly influenced by N1-H1ÁÁÁO2 interactions having the maximum E tot energy values; -51 and À27.2 kJ/mol for both experimental and predicted crystal structures. Next to N1-H1ÁÁÁO2 interaction, HÁÁÁH interactions in the experimental crystal structure and OÁÁÁH interactions in the predicted crystal structure contribute significantly over the crystal stability and the corresponding E tot energy values were À16.7 and À18.4 kJ/mol, respectively. The E tot energies calculated from crystal explorer are listed in Supplementary Table 1.

2-AAPA-NA complex
Molecular docking, the process of searching a ligand that is energetically and geometrically able to fit in the active site of a protein, is being widely used to discover novel drugs. In the present study, the optimized molecular structure of 2-AAPA was used for the molecular docking in the active site of NA protein. Molecular docking was carried out to understand the binding affinity, molecular conformation, and orientation of 2-AAPA molecules in the active site of NA protein.

Topological properties of electron density of hydrogen bonding interactions
The insight into the affinity of 2 AAPA conformers towards amino acid residue was given by the topological properties (  Figure 6a in which the red, white and blue surface denotes the hydrogen bond interaction, vdW separation and short contacts between the protein and the ligand. The sharp peaks in the fingerprint plot ( Figure  6b) confirms the strong hydrogen bon interactions between the protein and the ligand. The reduced density gradient [RDG] and its corresponding isosurface was obtained from Multiwfn (Lu & Chen, 2012) software. RDG is calculated to identify the NCI regions. The RDG scatter plot is shown in Figure 7a, which blue, green and red regions reveals the hydrogen bonding, van der Waals, and steric effect interactions between protein and ligand and the corresponding RDG surface is shown in Figure 7b.

Structural stability and flexibility of the docked complex
Simulation of molecular dynamics is an intuitive method for studying the conformation shifts of the NA-2-AAPA complex as a function of time in the presence of water solvent and is primarily necessary to obtain information on the stabilization of NA-2-AAPA interactions. Hence the docked complex is therefore subjected to simulations of molecular dynamics for 135 ns to examine the effect of the NA-2-AAPA interactions. Analyses such as RMSD, RMSF of individual residues,  hydrogen-bond occupancy, and secondary structure analyses have been analyzed in-depth to examine the interaction between the NA and 2-AAPA molecules.
To determine the effect of the 2-AAPA drug molecule on the conformational stability of the NA receptor protein, the 135 ns RMSD analysis corresponding to Ca, backbone, and all-atom of NA receptor protein was calculated for 135 ns concerning the starting structure of MD simulation and their corresponding results are presented in Figure 8a. There are lower values for the Ca and backbone RMSD than for the allatom RMSD. The NA-2-AAPA complex trajectory ranges from 0.350 to 2.175 Å, 0.403 to 2.224 Å, 0.435 to 2.474 Å, referring to the Ca, backbone, and all-atom RMSDs, respectively. During the simulation time of 0 to 13 ns, the Ca RMSD fluctuates from 0.350 to 1.309 Å. After that, the RMSD values rise with time changes and fluctuate between 1.381 and 2.175 Å during the simulation time from 13.1 and 36.8 ns. The NA-2-AAPA complex displays greater RMSD variations in the range of 1.916 to 2.175 Å over the 31.8 to 36.8 ns period, suggesting larger conformational shifts over this period. Beyond the 75 ns of simulation the RMSD variations gets saturated and fluctuates within the range of 1.597 to 2.031 Å and is held until the completion of the simulation. A similar trajectory pattern was found for the backbone and all-atom RMSD of the NA-2-AAPA complex. The average C alpha, backbone, and all-atom RMSDs of the NA-2-AAPA complex were found to be 1.664, 1.727, and 2.016 Å, suggesting the conformational stability of the complex, respectively.
To determine the structural flexibility of the NA protein, RMSF per protein residue was determined from the mean RMSF of its constituent residue and the corresponding findings are presented in Figure 8b. From the obtained results it is inferred that the residues VAL 83, ILE 84,ARG 107,LYS 111,LYS 143,ASN 146,THR 148,VAL 149,LYS 150,ARG 219,LYS 259,LYS 261,ASN 269,HIS 296,ASN 309,SER 342 .790 Å implying that flexibility of these residues whereas the remaining residues possess the smaller RMSF values in the range 0.390 to 1.980 Å which implies the rigidity of these residues. The rigidity of these residues is due to the formation of intramolecular hydrogen bonds. The superimposed view of the 2-AAPA-NA complex is presented in Figure 9. The structures were obtained from the initial (represented with pink color) and final structure (represented with blue color) of the NA-2-AAPA complex MD simulations and exhibit the conformational difference of ligands in the respective complex which allows us to visualize that how the conformation and the orientation of the ligands and proteins are altered during the MD simulation. However, the superimposed form of both docked and MD complexes do not exhibit many variations confirming that the 2-AAPA molecule binds more firmly with NA.

Occupancy of hydrogen bond analysis
The structural stability of the protein is characterized by hydrogen bond formation. Hence to characterize the dynamic behavior of the NA-2-AAPA complex, the occupancy of intra-and inter-molecular hydrogen bonds has been calculated and the obtained results are presented in Figures 10a,10b, and 11, respectively. In the case of intramolecular hydrogen-bond formation, the maximum and minimum occupancy of hydrogen bonds are found to be 329 and 269, respectively. The intra-molecular hydrogen-bond formation fluctuates between 288 and 329 during the first 20 ns. Beyond that the occupancy of hydrogen bonds decreases and fluctuates between 269 and 316, respectively. The inter-molecular hydrogen-bond formation between the NA-2-AAPA complex and solvent water molecules exhibits the maximum occupancy of hydrogen bonds in the range of 545-689. In this case, for the first 23 ns, the hydrogen-bond occupancy fluctuates between 545 and 645, and during the later stages, the occupancy of the hydrogen-bond increases and fluctuates between 649 and 689, respectively.

Variation of NA-2-AAPA interaction during the dynamics
In order to study the variation of the interactions between the NA protein and 2-AAPA ligand molecule during the simulation of 135 ns, the ligand-protein interactions were monitored at regular intervals of 5 ns from the starting structure (0 ns) to final structure (135 ns

Calculation of change in binding free energy using MM-GBSA method
In order to study the variation of binding free energy during the course of MD simulation, the Molecular Mechanics-Generalized Born model and Solvent Accessibility (MM-GBSA) energy calculations was carried out at regular intervals of 10 ns and also for the docked complex and final snapshot of MD simulation. The total binding free energy was contributed from the Coulomb energy, Covalent, H-bond, Lipo, p-p packing interaction, solvent generalized binding and van der Waals energy. The contribution of each term is presented in Table 8. The obtained data revealed that the binding free energy was found to be negative implying the stabilizing effect of the complex. Furthermore, the binding free energy decreases with increase in simulation time up to 80 ns and then it increases gradually up to À23.544 kcal/mol at 135 ns of MD simulation. This implies that during the final part of MD simulation the complex becomes more stable. Among all the interactions, the contribution of coulomb and van der Waals interactions to the total energy were found to be greater than the other interactions. The values of H-bond interaction of docked complex and the final snapshot of MD simulation were À0.025 kcal/mol, and À1.579 kcal/mol, respectively. It is obvious that after the simulation of 135 ns    the 2-AAPA drug molecule formed stable hydrogen bond with the residues of NA protein's binding pocket.

Conclusion
Neuraminidase is a potential drug target for the influenza virus. The ester-derived benzoic acid of 2-AAPA exhibits antiviral properties. Therefore, in this present study, X-ray crystallography, crystal structure prediction, Hirschfeld surface analysis, molecular docking, and MD simulation studies gave an insight into the molecular structure, strength of intermolecular interactions in the crystal phase, binding affinity, and stability of the molecule. The strength of the interaction between the molecules in the experimental and predicted crystal phase was highlighted from the analysis of HS and PIXEL energy calculations. Also, the topological analysis for closed-shell interactions between ligand and protein was carried out from Bader's Quantum theory of Atoms in Molecules. The stability of the 2-AAPA molecule in the active site of protein was evaluated from molecular simulations of 135 ns.