key: cord-0071246-zc08skoq authors: Jeyaram, R. A.; Radha, C. Anu title: Investigation on the Binding Properties of N1 Neuraminidase of H5N1 Influenza Virus in Complex with Fluorinated Sialic Acid Analog Compounds—a Study by Molecular Docking and Molecular Dynamics Simulations date: 2021-12-09 journal: Braz J Phys DOI: 10.1007/s13538-021-01009-z sha: b643430c2d008533064cc9c86782d48d8af695f7 doc_id: 71246 cord_uid: zc08skoq Worldwide, only two types of antiviral inhibitors (M2 ion channel protein inhibitor and Neuraminidase inhibitors) are approved to treat the influenza viral infection. But the mutation of amino acid sequence in the viral membrane proteins creates the viral resistance to existing antiviral drugs or inhibitors. So the corresponding antiviral drugs have to be reformulated to match these antigenic variations. Fluorination on the carbon–based molecule significantly enriches its biological properties. Hence this study is motivated to design the fluorinated sialic acid (SIA) analog inhibitors for the neuraminidase of H5N1 influenza A virus by substituting fluorine atom at different hydroxyls (O2, O4, O7, O8, and O9) of sialic acid. 100 ns molecular dynamics simulations are carried out for each protein–ligand complex system. NAMD pair interaction energy and MM–PBSA binding free energy calculations predict two possible binding modes for N1–SIA_F2, N1–SIA_F4, and N1–SIA_F7 complexes and single binding mode for N1–SIA_F8 and N1–SIA_F9 complexes. RMSD, RMSF, and hydrogen bonding analyses are used to understand the conformational flexibility and structural stability of each complex system. It has been concluded that the fluorinated sialic acid drug candidates SIA_F2 and SIA_F7 have better inhibiting potency against the N1 neuraminidase of H5N1 influenza virus. Sialic acid (SIA, Neu5Ac) is a nine carbon monosaccharide, which exists in more than 80 diversified forms due to the functional group modifications at C5 position and different hydroxyl oxygen's [1, 2] . Mostly SIA is present as the terminal residue of oligosaccharide chain of glyconjugate by making (2 → 3), (2 → 6) linkages with other glycans and (2 → 8), (2 → 9) linkages with other forms of sialic acids [3, 4] . SIA plays crucial roles in biological processes such as protein-carbohydrate interactions, recognition phenomena including cellular and molecular processes, communication, and signaling between the cells and molecules, cellular aggregation and development processes, tumor growth, and metastasis [5, 6] . Also SIA present on the host cell membrane acts as the receptor ligand for the host-pathogen (viruses, bacteria's, protozoa, and toxins) interactions. Most of the viral infections (Influenza virus, parainfluenza virus, murine hepatitis virus, noro virus, rota virus, DNA tumor viruses, and bovine corona virus) on the host cell can be initiated by the recognition of differently substituted acetyl group at different hydroxyl group of SIA [7, 8] . In the present work, we have given specific concern to a highly pathogenic virus, the H5N1 avian influenza virus. Influenza, generally referred to as flu, is an infectious disease caused by RNA viruses that infects both birds and mammals [9, 10] . Seasonally, this flu virus undergoes antigenic variations by altering their cell surface glycoproteins: hemagglutinins (HA) and neuraminidases (NA) causes the influenza epidemics [11, 12] . Globally, about 20% of the people are infected by these influenza epidemics in every year [13] . There are two types of antiviral inhibitors (M2 ion channel inhibitor and NA inhibitors) that are used to treat influenza viral infection, but the mutation in the viral structures causes the viral resistance to developed antiviral drugs or inhibitors [14] [15] [16] . So the corresponding antiviral drugs have to be reformulated to match these antigenic variations. Of the two important surface membrane proteins present in Influenza A viruses, HA is responsible for binding to the cell surface leading to the attachment and subsequent penetration by viruses into the target cell. NA is responsible for cleaving the terminal SIA moieties from the receptors to facilitate the elution of the progeny virions from the infected cell [17] [18] [19] . Hence, these two proteins can be the targeted for the design of inhibitors that can inhibit viral attachment and viral entry against influenza A virus. These two proteins recognize the SIA containing oligosaccharides present on the membrane of host cell receptors for their interaction which is crucial for viral entry into the cells. According to the mutagenic analysis the residues of both HA and NA binding sites are conserved for most of the influenza A strains. Hence, targeting HA and NA might inhibit the viral infection at the initial and post stages of influenza infections respectively. Simultaneous targeting will play a dual role in avoiding influenza infection. Substitution of fluorine atoms in biomolecules is widely used to enhance the metabolic properties of the drug candidate and it can dramatically alter the biological, electrochemical, and physicochemical properties of the biomolecule [20] . Due to the impact of the high electronegative charge, small molecular size and C-F bond stability; the fluorine atom significantly enhances the bioavailability, membrane and cell permeation, pKa, molecular lipophilicity, conformation of the molecule, and tissue distribution [21] [22] [23] [24] . In our earlier studies [25, 26] , we designed fluorinated SIA analog inhibitors for H5 HA of H5N1 influenza A virus by molecular dynamics (MD) simulation and we observed that fluorine substitution at O2, O7 and O9 hydroxyls of SIA increases the binding affinity towards the H5 protein. More recently, we have carried out molecular docking and MD simulation studies on the N1 NA of H5N1 influenza A virus complexed with zanamivir & sialic acid [27] . The catalytic center of NA enzyme makes stronger electro-statistical interaction to the SIA receptor and we observed that the SIA has much better inhibiting potency against N1 NA of H5N1 influenza A virus. In the continuation of our earlier work, the present investigation is focused on exploring the role of fluorination on different hydroxyl positions (O2, O4, O7, O8, and O9) of SIA complexed with N1 NA of H5N1 influenza A virus through molecular docking and MD simulation studies. In our previous investigation [27] , molecular docking was performed with sialic acid into the active site of N1 NA of H5N1 influenza A virus (pdb id: 3CKZ). In the present study, in-house developed Fortran program is used to substitute the fluorine atom into the different hydroxyl (O2, O4, O7, O8, and O9) positions of SIA leading to the formation of five different SIA analog molecules. The detailed drug likeliness, pharmacokinetics and toxicities of the each fluorinated SIA molecules are evaluated by the SWISS ADME webserver [43] . Autodock 4.2 software [44] is used to dock all the fluorinated SIA molecules into the active site residues of N1 of H5N1 influenza virus. Autodock tools (ADT) are employed to add the polar hydrogen atoms and Kollman charges to the N1 protein structure. The Gasteiger charges are assigned to the ligand atoms and the molecular grid size of 74 × 74 × 74 is set into the xyz direction with the grid spacing of 0.375 Å. The Lamarckian Genetic Algorithm (GA 4.2) is used for molecular docking of ligand molecule into the binding site of target protein and it gives the top ten binding score for the each protein-ligand complex structure. The best binding pose from each complex system is taken for the further MD simulation process. Amber 16 software package [28] is used to carry out the 100 ns MD simulation of N1fluorinated SIA complex structures. The input preparatory files are prepared through the antechamber and tleap modules. Generalized amber force field (gaff) and ff14SB force field parameters are used to fluorinated SIA ligands and N1 protein respectively. The 10 Å solvate waterbox of TIP3P water molecules (11,806 water molecules) are added from the solvent library and counter ions (5 Na + ions) are added to neutralize each complex structure. The steepest descent algorithm and the conjugate gradient algorithm are used in the energy minimization. The non-bonded cutoff interaction value is fixed at 8.0 and the temperature is fixed at 300 K with constant pressure. Finally the MD simulations of each complex structure are carried out through the particle mesh Ewald molecular dynamics (PMEMD) module [29] and the output trajectories are collected for every picosecond. The atomistic level interactions of the output trajectories are analyzed through Amber Tools17, visual molecular dynamics (VMD), nanoscale molecular dynamic simulation (NAMD), and with the in-house developed FORTRAN programs. Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) calculation [30] is used to calculate the binding free energy of the complex system by the formula ∆G bind = ∆G complex -(∆G N1 + ∆G fluorinatedSIA ). The NAMD pair interaction energy calculation [31] is computed between the binding site residues of N1 NA and the fluorinated SIA ligands. Total energy, root mean square deviation (RMSD), and root mean square fluctuation (RMSF) are calculated by the CPPTRAJ module of Ambertools 17 [32] . Chimera [33] , Molscript [34] , and Raster [35] are used to rendering the complex structures. In this study, the SWISS ADME web server (http:// www. swiss adme. ch/) is used to estimate the molecular physicochemical properties, lipophilicity, pharmacokinetics, druglikeness, leadlikeness, and synthetic accessibility of the ligand molecule. Lipinski rule of five is used to check the druglikeness (adsorption, distribution, metabolism, and elimination) of the lead molecule and it state that the (i) molecular weight less than 500 g/mol, (ii) calculated octanol-water partition coefficient, such as log P is less or than equal to 5, (iii) number of hydrogen bond acceptor ≤ 10, (iv) number of hydrogen bond donor ≤ 5 and (v) the molar refractivity must be 40-130 [45] . In our study, all the ligand molecules passes the Lipinski rules with the violation of total number of hydrogen bond donor is greater than 5. Details of the ADME evaluations are given in Table 1 . It has been observed that all the fluorinated SIA molecules can be used a drug candidate to inhibit the target protein structure. The fluorinated SIA molecules are successfully docked into the active site of N1 NA of H5N1 Influenza virus by Autodock 4.2 software. The docked structure of each complex system is shown in the Fig. 1 and the estimated free energy for all the complex systems are given in Table 2 . The autodock program uses the formula [K i = exp(DeltaG/(R*T))] to convert the binding energy into the inhibition constant (K i ). This may be calculated further experiments in future. Both binding free energy and K i values of all the complex systems are very close to the docking of N1-SIA complex (binding free energy = -6.88 kcal/mol and K i = 31.02 μM) reported in our previous study [27] . The docking study predicted that the active site residues R118, E119, D151, R152, E276, E277, R292, N294, Y347, R371, and Y406 of N1 protein bind with the fluorinated SIA molecules which is similar to the binding site residues involved as reported in our earlier study [27] . Next, all the fluorinated SIA molecules complexed with N1 NA of H5N1 influenza viruses are computationally simulated by Amber16 software package and the output trajectories are thoroughly analyzed by the analytical tools. The binding stability and conformational flexibility of each complex system is analyzed by the RMSD, RMSF and the direct and water mediated hydrogen bonding interactions. NAMD pair interaction energy calculation, MM-PBSA binding free energy calculation and total energy calculations are used to identify the binding efficacy of each protein-ligand complex system. The NAMD pair interaction energy of all the complex systems are computed between the binding site residues (R118, E119, V149, D151, R152, R156, W178, N221, E227, G244, S246, N247, E276, E277, R292, N294, N325, A346, Y347, G348, N369, R371, Y406, R430, P431 and T437) of N1 protein and the fluorinated SIA ligands are given in Fig. 2 . NAMD predicts two possible binding modes (BM1 and BM2) for the N1-SIA_F2 (BM1 = 1-37 ns; BM2 = 38-100 ns), N1-SIA_ F4 (BM1 = 1-39 ns; BM2 = 40-100 ns) and N1-SIA_F7 It is observed that the BM1 of N1-SIA_F2 and N1-SIA_ F7 complexes have the highest binding efficacy towards N1 protein. The MM-PBSA binding free energy is calculated for all the complexes to validate the outcome of pair interaction energy. This also agrees well with the better inhibiting potency of N1-SIA_F2 and N1-SIA_F7 complexes in BM1 as predicted by NAMD. The binding energy difference between N1-SIA_F2 and N1-SIA_F7 is 13 kcal/mol in NAMD and 1 kcal/mol in MM-PBSA. The order of inhibiting potency of fluorinated SIA molecules against the N1 of H5N1 influenza A virus is SIA_F2 ≈ SIA_F7 > SIA_F9 ≈ SIA_F8 > SIA_F4 and the details of NAMD pair interaction energy and MM-PBSA binding free energies are given in Table 2 . In our earlier studies, we reported that the magnitude of energy difference between NAMD and MM-PBSA is due to the impossibility of consideration water mediated interaction in NAMD [25, 26] . Total energy (E tot ) is the sum of energies from E bond + E angle + E dihedral + E vdw14 + E elec14 + E vdw + E elec and it is calculated for each complex system by the CPPTRAJ module of AmberTools 17. The average total energy for each complex system is given in Table 2 , which agrees with the trend of binding efficacy predicted by the NAMD pair interaction energy and MM-PBSA binding free energy calculations. Visual molecular dynamics (VMD) software is used to analyze the hydrogen bonding and hydrophobic interactions between each protein-ligand complex system [36] . Hydrogen bonding interactions observed between the donor and acceptor atoms with the distance cutoff 3.5 Ǻ and angle cutoff 30°. From the hydrogen bonding analysis, it is observed that the active site of N1 NA keeps the SIA in 2 C 5 chair conformation and the schematic representation is given in Fig. 3 . Total number of hydrogen bonds in BM1 and BM2 and the details of interaction patterns in BM1 are provided in Tables 2 and 3 , respectively. The BM1 of N1-SIA_F2 complex makes 10 direct hydrogen bonds and 4 water mediated hydrogen bonds. O1 carboxyl group of SIA_F2 form the direct hydrogen bond with the side chain nitrogen of R292, R371 and side chain hydroxyl oxygen of Y347. The O4 ring oxygen and O10 acetamide oxygen interact with side chain nitrogen of R118. The acetamide nitrogen N5 makes direct hydrogen bond with the side chain hydroxyl of Y406. The glycerol side chain oxygen O8 interacts with the side chain of R292 and Y347. Water-mediated hydrogen bonding interaction is observed between the glycerol side chain oxygen O7 and side chain of nitrogen of R152. The glycerol side chain O8 make a water-mediated hydrogen bond with side chain of E277 and N294. At last, the acetamide oxygen O10 and side chain of R156 interacted through water mediated hydrogen bonds. All the water mediated hydrogen bonds having the average residual time of ≥ 6 ps. BM1 of N1-SIA_F4 complex form 5 direct hydrogen bonds (S165, R292, and N294) and 9 water mediated hydrogen bonds (D151, R152, E276, E277, R371 and G348) and BM1 of N1-SIA_F7 makes 8 direct (E277, R292, N294) and 8 water-mediated hydrogen bonds (D151, R152, R156, W178, and E277) to stabilize the complex structure. It is observed that the N1-SIA_F8 and N1-SIA_F9 complexes have only one binding mode (BM1). The N1-SIA_F8 complex forms 6 direct and 4 water mediated hydrogen bonds, after 90 ns the ligand moves around 70 Å away from the active site of N1 protein. In N1-SIA_F9 complex, 8 direct and 5 watermediated hydrogen bonds are observed to stabilize the complex structure. In the BM2 of N1-SIA_F2 (38-100 ns) complex only 5 direct hydrogen bonds (V149, R430, P431, and T437) and 2 water mediated hydrogen bonds (R118) are observed. Fluorine substitution on O2 (named as F2) directly interacts with side chain nitrogen of R430. The O4 ring oxygen makes direct hydrogen bonds with the backbone oxygen of R430, backbone nitrogen of P431, and side chain hydroxyl group of T437. The acetamide oxygen O10 forms water mediated hydrogen bonds with backbone oxygen of V149 and side chain nitrogen of R118. The N1-SIA_F4 complex in BM2 (40-100 ns) also forms 5 direct hydrogen bonds (N221, G244, S246, N247, and N294) and 2 water mediated hydrogen bonds (Y347). The BM2 of N1-SIA_F7 complex (44-49 and 64-100 ns) make 5 direct hydrogen bonds (R152, S246) and 4 water-mediated hydrogen bonds (R152, N221, and Y347). The structural stability of the binding modes is mainly due to the hydrogen bonds and also due to hydrophobic interactions between ligand-protein components. For BM1 of N1-SIA_F2, R371HH21/22, and Y406HD1 interact with H3 of SIA_F2, H4 makes the hydrophobic contact with R292HH22 and Y406HD1&HE1. H5A interact with Y406HE1/HE2 and H6 interact with R292HH22. There is no hydrophobic interaction observed in glycerol side chain hydrogen's of SIA_F2. The N-acetamide hydrogen's H11, H12, and H13 forms 2 to 4 hydrophobic interactions with E119HG2 and HG3 and E227HG2 and HG3. In BM2 of N1-SIA_F2, H3 forms the hydrophobic contacts (4 to 9 bonds) with V149, R430, P431, and T439. H5 makes 2 or 3 contacts with V149/I427/R430. The glycerol side chain hydrogen H7 interact with V149/P431, H8 interact with V149/R430, and H9 interact with V149/I47, P431, and T439. Finally, N-acetamide hydrogen's (H11, H12, and H13) interact with R118 and T439 or (V149, R430, P431, and T439). The BM1 of N1-SIA_F4 makes the maximum of 15 hydrophobic interactions with the N1 protein (I2222, R224, S246, R292, N294, and Y347). On the other hand, BM2 of N1-SIA_F4 forms the maximum of 25 hydrophobic interactions with the amino acid residues (N221, I222, R224, P245, S246, N247, N294, and Y347). It is observed that more number of hydrophobic clusters are observed in BM2 of N1-SIA_F2 and N1-SIA_F4. So that the ligand remains in the same binding mode to maintain the complex form and cannot return back to the BM1. This might be the reason for the longer residential time of BM2 of these complexes when comparing to its BM1 binding mode in NAMD calculation. BM1 of N1-SIA_F7 makes the hydrophobic interactions with E119, D151, R152, E227, S246, R292, N294, Y347, R371 and Y406. The BM2 of N1-SIA_F7 interacts with R152, N221, I222, P245, S246, N247, R292, and Y347. The N1-SIA_F8 makes the hydrophobic interactions with R118, D151, R152, E227, S246, R292, N294, Y347, and Y406. The N1-SIA_F9 interacts with R292, P326, Y347, R371, and P431. In protein-ligand complex system, after the MD simulation, the RMSD and RMSF are used measure the binding stability and conformational flexibility of the protein structure. RMSD is a measure of average deviations from the initial structure to its final conformation. In this study, RMSD is measured for the protein backbone Cα atom with respect to the simulation time by using CPPTRAJ module of Amber16. In all the complexes, the protein structure is deviated from 0.5 Å to 2.3 Å. RMSF is used to measure the conformational flexibility of the protein (backbone Cα atom) structure based on the residual fluctuation. It is observed that the protein structures are fluctuates up to 6 Å and the large RMSF fluctuations are observed in the region of 142-150, 340-343, 430-434, and at the end of the amino acid residues, which are not corresponds to the binding site residues and these fluctuations are due to flexibility of these loop regions. The RMSD and RMSF values are given in Fig. 4 . The same trend of deviations and fluctuations are noticed in our earlier study [27] and the similar trend of measurements was reported in the inhibitor design for N1 NA of H5N1 and H1N1 influenza A viruses [37] [38] [39] . In the extension of our previous study, we docked the fluorinated SIA molecules into the active site of N1 NA of H5N1 influenza A virus and all the N1 protein-fluorinated SIA ligand complexes are simulated by Amber16. Interestingly, the N1 active site accommodate all the fluorinated SIA molecules in its solution state 2 C 5 conformation like N1-SIA complex reported in our previous study [27] . RMSD and RMSF plots show the conformational flexibility and binding stability of each complex system. The NAMD pair interaction energy, MM-PBSA binding free energy and total energy calculations agree well with one another and it reveals that the better inhibiting potency of N1-SIA_F2 and N1-SIA_F7 complexes in BM1. In our earlier study, we report that the SIA_F2 and SIA_F7 drug models have the better inhibiting potency against H5 HA of H5N1 influenza virus [25] . [27, [40] [41] [42] . In our previous investigation, we reported that N1 neuraminidase adopt wide variety of conformations due to the flexibility of 150 loop (the amino acid residues 147-152) forms a large cavity adjacent to the active site residues and it allows the second SIA binding site residues (S370, S372, I400, T401, W403, K432) interaction along with the active site of N1 NA of H5N1 influenza virus [27] . But the fluorinated SIA molecules in BM1 avoid the presence of second SIA binding site interaction; this may be the reason for the weaker inhibiting potency of the N1-SIA_F4, N1-SIA_F8 and N1-SIA_F9 complexes. From the Table 1 , it is clear that the BM2 of N1-SIA_F2, N1-SIA_F4, and N1-SIA_F7 complexes have 5 direct and 2 to 4 water-mediated hydrogen bonds which reflect the reduced energy contribution in NAMD and MM-PBSA calculation. All the analyses conforms that the higher binding efficacy of N1-SIA_F2 and N1-SIA_F7 complexes in its BM1 and the order of inhibiting potency against N1 NA is SIA_F2 ≈ SIA ≈ SIA_F7 > SIA_F9 ≈ SIA_F8 > SIA_F4. In the continuation of our earlier studies (fluorinated SIA molecules inhibitor for H5 HA of H5N1 influenza A virus), we designed the fluorinated SIA analog inhibitors for N1 NA of H5N1 influenza A virus by using molecular docking and MD simulations. RMSD, RMSF, and hydrogen bonding analyses are used to understand the binding stability and conformational flexibility of each complex system. NAMD pair interaction and MM-PBSA calculations predict two possible binding modes for N1-SIA_F2, N1-SIA_F4, and N1-SIA_F7 and single binding modes for N1-SIA_F8 and N1-SIA_F9. We suggest that one can design and use SIA_F2 and SIA_F7 inhibitors or antiviral drugs to inhibit the viral activity of H5 HA and N1 NA of H5N1 influenza A virus. In the future, it is possible to design fluorine based inhibitors for different viral diseases. Sialic acid and biology of life: an introduction', sialic acids and sialoglycoconjugates in the Biology of Life, Health and Disease Molecular virology of human pathogenic viruses Influenza other Respir. Viruses Raster: Geographic analysis and modeling with raster data The authors declare that there is no conflict of interest in this work.