Structure-based virtual screening of influenza virus RNA polymerase inhibitors from natural compounds: Molecular dynamics simulation and MM-GBSA calculation

Computational Biology and Chemistry 85 (2020) 107241

Contents lists available at ScienceDirect

Computational Biology and Chemistry journal homepage: www.elsevier.com/locate/cbac

Research Article

Structure-based virtual screening of influenza virus RNA polymerase inhibitors from natural compounds: Molecular dynamics simulation and MM-GBSA calculation

T

Zhe Jina, Ying Wanga, Xiao-Fei Yua, Qi-Qi Tana, Shi-Shao Lianga, Tai Lia, Hong Zhangb, Pang-Chui Shawc, Jian Wanga,*, Chun Hua,* a b c

Key Laboratory of Structure-Based Drug Design & Discovery, Ministry of Education, Shenyang Pharmaceutical University, Shenyang, 110016, Liaoning, China School of Life Sciences and Biopharmaceuticals, Shenyang Pharmaceutical University, Shenyang, 110016, Liaoning, China School of Life Sciences, The Chinese University of Hong Kong, Hong Kong SAR, China

ARTICLE INFO

ABSTRACT

Keywords: Influenza virus Protein-protein interaction inhibitor Virtual screening Drug likeness Molecular dynamics simulations

The resistances of matrix protein 2 (M2) protein inhibitors and neuraminidase inhibitors for influenza virus have attracted much attention and there is an urgent need for new drug. The antiviral drugs that selectively act on RNA polymerase are less prone to resistance and possess fewer side effects on the patient. Therefore, there is increased interest in screening compounds that can inhibit influenza virus RNA polymerase. Three natural compounds were found by using molecular docking-based virtual screening, which could bind tightly within the polymerase acidic protein-polymerase basic protein 1 (PA-PB1) subunit of influenza virus polymerase. Firstly, their drug likeness properties were evaluated, which showed that the hepatotoxicity values of all the three compounds indicating they had less or no hepatotoxicity, and did not have the plasma protein biding (PPB) ability, the three compounds needed to be modified in some aspects, like bulky molecular size. The stability of the complexes of PA-hits was validated through molecular dynamics (MD) simulation, revealing compound 2 could form more stable complex with PA subunit. The torsional conformations of each rotatable bond of the ligands in PA subunit were also monitored, to investigate variation in the ligand properties during the simulation, compound 3 had fewer rotatable bonds, indicating that the molecule had stronger rigidity. The bar charts of protein–ligand contacts and contacts over the course of trajectory showed that four key residues (Glu623, Lys643, Asn703 and Trp706) of PA subunit that participated in hydrogen-bond, water bridge and hydrophobic interactions with the hit compounds. Finally, the binding free energy and contributed energies were calculated by using MM-GBSA method. Out of the three compounds, compound 1 showed the lowest total binding free energy. Among all the interactions, the contribution of the covalent binding and the van der Waals energy were more than other items, compound 1 formed more stable hydrogen bonds with the residues of PA subunit binding pocket. This study smoothed the path for the development of novel lead compounds with improved binding properties, high drug likeness, and low toxicity to humans for the treatment of influenza, which provided a good basis for further research on novel and effective influenza virus PA-PB1 interaction inhibitors.

Abbreviations: M2 protein, Matrix protein 2; RNA, ribonucleic acid; PA, polymerase acidic protein; PB1, polymerase basic protein 1; ADMET, absorption distribution, metabolism, excretion, and toxicity; HBA, hydrogen-bond acceptor; HBD, hydrogen-bond donor; TPSA, topological polar surface area; IAV, influenza A virus; IBV, influenza B virus; ICV, influenza C virus; PSE, ethanolic extract of Peperomia sui; NP, nucleoprotein; TGBG, 13,4,6-tetra-O-galloyl-β-D-glucopyranoside; PB2, polymerase basic protein 2; FDA, food and drug administration U.S; mRNA, messenger ribonucleic acid; RCSB, research collaboratory for structural bioinformatics U.S; PDB, protein data bank; PAC, C-terminal region of PA; PB1N, N-terminal region of PB1; SP, standard precision; MW, molecular weight; Fraction Csp3, fraction of carbons in the sp3 hybridization; GI, gastrointestinal; BBB, blood brain barrier; P-gp, P-glycoprotein; PSA, polar surface area; PPB, plasma protein binding; MD, Molecular dynamics; RMSD, root mean square deviation; rGyr, radius of gyration; intraHB, intramolecular hydrogen bond; MolSA, molecular surface area; SASA, solvent accessible surface area; VS, virtual screening; MM-GBSA, molecular mechanics / generalized born surface area ⁎ Corresponding authors. E-mail addresses: [email protected] (J. Wang), [email protected] (C. Hu). https://doi.org/10.1016/j.compbiolchem.2020.107241 Received 13 October 2019; Received in revised form 22 February 2020; Accepted 26 February 2020 Available online 26 February 2020 1476-9271/ © 2020 Published by Elsevier Ltd.

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al.

1. Introduction Influenza, commonly known as the flu, is an infectious disease caused by the influenza virus. Two of the four types of influenza viruses affect humans: Type A (influenza A virus, IAV) and Type B (influenza B virus, IBV). The structures of IAV and IBV are very similar, the virus particle is elliptical and 80–120 nanometers in diameter (Sugita et al., 2011). The structure of influenza C virus (ICV) is filamentous, which is different from the shape of IAV and IBV. However, despite there are varieties in shapes, the viral particles of all influenza viruses are quite similar in composition (Bouvier and Palese, 2008). Hemagglutinin (HA) and neuraminidase (NA) are the two large glycoproteins on the outside of the flu viral particles (Suzuki, 2005). Furthermore, they are antigens to which antibodies can be raised. Influenza A viruses are classified into different subtypes based on the antibody responses to HA and NA. These different types of HA and NA form the basis of the H and N subtypes distinctions in, for instance, H1N1 and H3N2 (Hilleman, 2002). There are 18 H and 11 N subtypes, but H1N1 and H3N2 subtypes are main cause of disease in humans (Tong et al., 2013; Wu et al., 2014). IBV consists of B/Yamagata lineage and B/Victoria lineage also can cause the human disease (van de Sandt et al., 2015). Currently, three classes of antiviral influenza virus drugs are widely used in clinic. They are the M2 protein inhibitors (amantadine derivatives), the neuraminidase inhibitors (oseltamivir, zanamivir, peramivir and laninamivir) and the cap-dependent endonuclease inhibitor (baloxavir marboxil). The M2 protein inhibitors are not effective against IBV because of the virus lack the protein (Kaplan, 2010). During the 2005–2006 influenza season, the global incidence of resistance to adamantanes among A(H1N1) isolates was 15.5 %, accounting for 66 of 424. Among them, in China, the incidence of resistance was 71.7 % (33/46); in Asia, the incidence was 22.6 % (45/199). However, during the same influenza season, the global incidence of resistance to adamantanes among A(H3N2) isolates was 90.6 %, accounting for 1059 of 1169. Among them, in China, the incidence of resistance was 100 % (8/8); in Asia, the resistance was 78.5 % (113/144) (Deyde et al., 2007). In Foshan district of China, in 2008, the incidence of resistance to adamantane among A(H1N1) isolates was 92.2 % (59/64), and in 2009, the value was 31.3 % (26/83). Oseltamivir resistance was only 4.4 % (6/137) in 2008, and in 2009, the value became 94.4 % (101/107) (Zhou et al., 2011). Since then, the resistance of M2 protein inhibitors and neuraminidase inhibitors have attracted more attention. Many natural products showed inhibitory activity on influenza virus. The ethanolic extract of Peperomia sui (PSE) significantly increased the viability of cells infected by the H6N1 virus. PSE also suppressed the synthesis of viral nucleoprotein (NP), and inhibited the growth of the virus in DF-1 cells. Further, PSE was found to inhibit the neuraminidase activity of H6N1 virus (Yang et al., 2014). Dendrobine, a major component of Dendrobium nobile, showed antiviral activity against influenza A viruses, including A/FM-1/1/47 (H1N1), A/Puerto Rico/8/34 H274Y (H1N1), and A/Aichi/2/68 (H3N2) with IC50 values of 3.39 ± 0.32, 2.16 ± 0.91, 5.32 ± 1.68 μg/mL, respectively. Dendrobine could bind to the highly conserved region of viral NP, subsequently restraining nuclear export of viral NP and its oligomerization (Li et al., 2017b). 1,3,4,6-Tetra-O-galloyl-β-D-glucopyranoside (TGBG, Fig. 1 a) from Euphorbia humifusa Willd showed broad-spectra anti-influenza activity against two seasonal influenza A strains, A/California/ 07/2009 (H1N1) and A/Perth/16/2009 (H3N2), and seasonal influenza B strain B/Florida/04/2006. Immunofluorescence assay demonstrated that TGBG significantly inhibited nuclear export of influenza NP during the early stages of infection, causing NP to accumulate in the nucleus. In addition, influenza-induced activation of the Akt signaling pathway was suppressed by TGBG in a dose-dependent manner (Chang et al., 2016).

Fig. 1. The natural compounds that showed the inhibitory activity on influenza virus, 1,3,4,6-tetra-O-galloyl-β-D-glucopyranoside, TGBG(a), flutamide (b) and aromatic analogues of flutimide (c).

The natural product flutimide (Fig. 1 b) was isolated from Delitschia confertaspora. It is the first natural product shown to inhibit the capdependent endonuclease of influenza virus, and the IC50 value of endonuclease inhibition was 3 μM (Hensens et al., 1995). Some aromatic analogues of flutimide were found to exhibit good inhibitory activity on influenza virus, and the most potent compound with p-methoxybenzylidene group at C-5 of 3H-pyrazine-2,6-dione (Fig. 1 c) showed IC50 values of 0.8 μM (Singh and Tomassini, 2001). With the growing problem of resistance, there is an urgent need to focus on new other targets. RNA polymerase is an important target against influenza virus. The polymerase consists of polymerase acidic protein (PA), polymerase basic protein 1 (PB1) and polymerase basic protein 2 (PB2) subunits. The cap-dependent endonuclease inhibitor xofluza (baloxavir marboxil) was approved by U.S. FDA in Oct. 2018 (U.S. FDA, 2018). The drug is the first new antiviral flu treatment with a novel mechanism of action approved by the FDA in nearly 20 years. The polymerase may synthesize viral messenger RNAs (mRNAs) using short capped primers derived from cellular transcripts by a unique ‘capsnatching’ mechanism. The PB2 subunit binds the 5′ cap of host premRNAs, which are subsequently cleaved after 10–13 nucleotides by the viral endonuclease. The biochemical and structural studies show that 209 residues at the N-terminal of the PA subunit contain the endonuclease active site, and conclusively the anti-influenza activity is inhibited by the inhibitor of the influenza cap-dependent endonuclease (Dias et al., 2009). In 2008, there were two X-ray crystal structures on IAV RNA polymerase PA-PB1 complex reported. The RCSB PDB codes were 2 ZN L (Obayashi et al., 2008, Fig. 2 a) and 3CM8 (He et al., 2008, Fig. 2 b). 2 ZN L derived from the complex of PA (C-terminal region of PA, PAC, residues 239–716) and PB1 (N-terminal region of PB1, PB1N, residues 1–81), Flu A virus strain A/Puerto Rico/8/1934 H1N1 at 2.3 Å resolution. 3CM8 derived from the complex of PA (PAC, residues 257–716) and PB1 (PB1N, residues 1–25), avian Flu A virus strain A/ Hong Kong/156/1997 H5N1at 2.899 Å resolution. The PB1N and PAC of the RNA polymerase interaction is essential viral transcription and replication. The antiviral drugs that selectively act on RNA polymerase are less prone to resistance and have less adverse effect on human. The IAV RNA polymerase PA-PB1 complex, with the PDB codes 2 ZN L (a) and 3CM8 (b) are shown in the supplementary materials. The previous research on influenza virus suggests that, the interactions between PA and PB1 are inhibited, the function of influenza virus RNA polymerase is affected, and the activity of influenza virus is inhibited (Ng et al., 2009; Massari et al., 2013; Pagano et al., 2013; Liu et al., 2015, 2017a). To find PA-PB1 interaction inhibitor from natural source, structure-based virtual screening was performed based on the PDB file 2ZNL derived from IAV H1N1 subtype. 2. Materials and methods 2.1. Structure-based virtual screening Molecular docking was a crucial technique for exploring ligandprotein interactions (included such as hydrogen-bond interaction, π-π 2

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al.

interaction, hydrophobic interaction) and calculating the corresponding binding free energy (Wang et al., 2016). The standard precision (SP) Glide algorithm (Friesner et al., 2006) was used for the molecular docking in our previously published work (Wang et al., 2018a), and the algorithm had been proven to be accurate and reliable. The purpose of docking was to explore the interactions between PA subunit of the RNA polymerase and small molecules ligand. The crystal structure of IAV RNA polymerase PA-PB1 complex (2ZNL.pdb) was retrieved from RCSB Protein Date Bank (RCSB PDB, http://www.rcsb.org/) (Obayashi et al., 2008). The downloaded IAV RNA polymerase PA-PB1 complex protein file (2ZNL.pdb) was prepared by Schrödinger’s protein preparation wizard, which removed PB1 subunit and water molecules from the protein, built the missing side chain atoms, corrected the structure of the protein and added hydrogen atoms before docking (Shivakumar et al., 2010). Afterwards, the OPLS_2005 force field was used to optimize protein energy and eliminate steric hindrance (Itteboina et al., 2017). The position of PB1 subunit was defined to generate grid for subsequent molecular docking work (Wang et al., 2018a). The grid box of 2ZNL cavity was generated, and then the molecules were flexibly docked into the prepared protein using the XP precision in Glide. The ligand docking panel implemented in 2015 Schrödinger software package was employed to dock the ligands to IAV RNA polymerase PA subunit with three modes of Glide docking (Aier et al., 2016). The parameters in the docking process were set as defaults. The protonation states of residues were determined at physiological pH 7.4, which was performed via the protein preparation wizard within Schrödinger software package. A total of 70,000 compounds available for virtual screening were obtained from some natural compounds databases. Three steps of virtual screening, including HTVS, SP and XP involved Glide score for representing binding affinity and ranking ligands. HTVS mode of docking was performed to conformers of 70, 000 compounds and the top 10 % ranking ligands were further selected for subsequent SP docking, and finally the top 10 % were re-docked with subsequent Glide XP docking. The screened three compounds with high docking scores were obtained through this rigorous docking based in-silico filtering. Finally, the docking scores, hydrogen bond and VDW interactions were visualized in docking report which was used for further analysis. 2.2. Lipinski’s rule and ADMET prediction In this work, the physicochemical properties and different pharmacokinetics parameters of the screened compounds were calculated using Swiss ADME prediction (http://www.swissadme.ch/) and Toxicity prediction module of Discovery Studio v3.0 (Accelrys, USA) to calculate Lipinski’s rule-of-five drug-likeness properties for oral bioavailability and ADMET properties to determine the modification direction of the compounds. Firstly, Swiss ADME online server (http://www.swissadme.ch/) was used to calculate the ADME properties of the investigated compounds. The following data for compounds 1, 2 and 3 are shown in the bioavailability radar charts: lipophilicity (LIPO), molecular weight (SIZE), polarity (POLAR), solubility (INSOLU), saturation (INSATU) and flexibility (FLEX). The following data were shown in the data table: molecular weight (MW), fraction of carbons in the sp3 hybridization (Fraction Csp3), number of rotatable bonds, number of hydrogen-bond acceptors (HBA), number of hydrogen-bond donors (HBD), TPSA, Log Po/w (XLOGP3), log S (ESOL), solubility, gastrointestinal (GI) absorption, blood brain barrier (BBB) permeation, P-glycoprotein (P-gp) substrate, cytochrome P450 enzyme inhibitory capacity (1A2, 2C9, 2C19, 2D6, 3A4 subtypes), skin permeation (log Kp). Secondly, ADMET descriptors algorithm and Toxicity prediction (extensible) module of Discovery Studio v3.0 (Accelrys, USA) were used to calculate the physicochemical properties, drug-likeness and different pharmacokinetics parameters of the compounds (Liu et al., 2019). The results of Lipinski’s rule calculation for compounds 1, 2 and 3,

Fig. 2. Two-dimensional (2D) diagrams of ligand-protein interactions in the active sites of the PA subunit of IAV RNA polymerase (2 ZN L.pdb), compound 1 (a), 2 (b) and 3 (c). The hydrogen-bond interactions with residues are represented by a purple dashed arrow directed towards the electron donor. The ππ interactions are represented by a green or a red line.

3

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al.

include MW, number of HBA, number of HBD, polar surface area (PSA), QP log S were obtained from Discovery Studio prediction. Cytochrome P450 2D6 (CYP2D6) enzyme inhibitory capacity, hepatotoxicity, plasma protein binding (PPB) ability of compounds 1, 2 and 3 were calculated by ADME prediction.

Table 1 The structures and virtual screening results of molecules. Number

2.3. Molecular dynamics simulations Molecular dynamics (MD) was performed on Desmond v4.3 suite, running on Red Hat 6.10 Linux operating system. Noose-Hover chain thermostat at 300 K, Martyna-Tobias-Klein barostat at 1.01325 bar, isotropic coupling, Coulombic cutoff at 0.9 nm. The rest of the parameters were default. First, the protein-ligand complexes (crystal structure of PA subunit bound with three lead compounds obtained by docking) were prepared in protein preparation wizard of Desmond module with the explicit solvent model and SPC and orthorhombic box shape. The box size was set at distance of 10 Å from the outermost atom of the protein-ligand complex. The sodium chloride with the approximately physiological concentration of 0.15 M was placed to the simulation box to set the ionic strength. After the minimization jobs were performed to relax the system into a local energy minimization, this model system was submitted to MD simulation steps using Desmond v4.3 package with OPLS_2005 force field. The specific simulations steps were operated similar with the previous work (Wang et al., 2018b; Chen et al., 2018; Liu et al., 2019; Wang et al., 2019). In Maestro, the quality of MD simulations was assessed by the Simulation Quality Analysis tools and analyzed by the Simulation Event Analysis tool, and ligand-receptor interactions were identified using the Simulation Interaction Diagram tool.

Structure

Interaction with 2 ZN L

MM/ GBSA (kcal/ mol)

Dock score (kcal/ mol)

1

Phe411 Gly622 Glu623 Thr639 Lys643 Asn703 Trp706 Ser709

−77.00

−10.752

2

Asn412 Glu623 Glu630 Lys643 Asn703

−58.18

−10.574

3

Phe411 Cys415 Glu623 Lys643 Asn703

−68.64

−10.452

2.4. Binding free energy calculation The docked conformations were energy-minimized by the Prime (version 3.5) module, and then Molecular Mechanics Generalized Born Surface Area (MM-GBSA) was applied to calculate the binding free energy of the complexes. The OPLS_2005 force field (Shivakumar et al., 2010), GBSA continuous solvent model and the rotamer search algorithms were used throughout the calculation.

having the similar binding affinity with the active sites of the PA subunit, were all regarded as representatives for detailed description. Fig. 3 shows that compounds 1, 2, and 3 basically occupied the same binding pocket of the PA subunit of IAV RNA polymerase. There are two adjacent amino acid residues, Lys643 and Asn703, located at one side of the binding pocket; and residue Glu623 on the other side. The molecular docking results showed compounds 1, 2, and 3 generated nearly identical hydrogen bond interactions with the three important amino acid residues. Meanwhile, there were some differences among them. Compounds 1 and 3 generated similar interactions with amino acid residues around the two compounds. There were π-π interactions between the compounds and residue Phe411. The two compounds generated hydrophobic interactions with the amino acid residues around these compounds, Gln408-Phe411-Asn412-Cys415, and Phe707-Ser709-Phe710, respectively. There were no similar interactions found between compound 2 and the above-mentioned residues.

GGB = EMM + Gsolv + GSA In the equation, the electrostatic solvation energy (△GGB) was calculated by using the general borne (GB) Models; △EMM was the difference between the minimized energies of ligand-protein complexes and the total energies of protein and ligand in free form; △Gsolv was the difference in the GSA solvation energies of the ligand-receptor complex and the sum of the solvation energies of receptor and the ligand in the unbound state; △GSA is the difference in the surface area energies for the protein and the ligand. 3. Results and discussion 3.1. Virtual screening hits Based on the docked poses within the active sites of the PA subunit of IAV RNA polymerase (2 ZN L.pdb), top three molecules with lower docking score were identified. Molecular docking and binding energies were predicted through CDOCKER within Accelrys Discovery Studio 2.5.5 package. The coordinate of the PA–compound complex (PDB code: 2ZNL) (Obayashi, et al.) was employed. The results are listed in Table 1. For each compound, the docked conformation which had the lowest docking score was selected to analyze the binding pattern. All these compounds had 1 or 2 interactions with any of the key three amino acids (Glu623, Lys643 and Asn703). Out of the top 3 molecules, compound 1 (Glide score: -10.752 kcal/mol), compound 2 (Glide score: -10.574 kcal/mol), and compound 3 (Glide score: -10.452 kcal/mol),

3.2. Lipinski’s rule and ADMET prediction Two methods were used to calculate Lipinski’s rule-of-five druglikeness properties for oral bioavailability and ADMET properties and to determine the direction of modification. Firstly, Swiss ADME prediction (http://www.swissadme.ch/) was used to calculate Lipinski’s rule-of-five drug-likeness properties. In Figure 11, The pink area represents the optimal range for each property. LIPO stands for lipophilicity, XLOGP3 value should be between -0.7 and +5.0; SIZE, molecular weight should be between 150 and 500 4

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al.

Fig. 3. The bioavailability radar charts of compounds 1 (a), 2 (b) and 3 (c) obtained from Swiss ADME prediction.

Daltons; POLAR stands for polarity, topological polar surface area (TPSA) should be between 20 and 130 Å2; INSOLU stands for solubility, log S should be not higher than 6; INSATU stands for saturation: fraction of carbons in the sp3 hybridization should be not less than 0.25; and FLEX stands for flexibility, there are should be no more than 9 rotatable bonds (Daina et al., 2017). The results of Lipinski’s rule calculation for compounds 1, 2 and 3, include molecular weight (MW), fraction of carbons in the sp3 hybridization (Fraction Csp3), number of rotatable bonds, number of Hbond acceptors (HBA) should be not less than 10, number of H-bond donors (HBD) should be not less than 5, TPSA, Log Po/w (XLOGP3), log S (ESOL) that obtained from Swiss ADME prediction are listed in Table 2. The ADME prediction results for compounds 1, 2 and 3, include solubility, gastrointestinal (GI) absorption, blood brain barrier (BBB) permeation, P-glycoprotein (P-gp) substrate, cytochrome P450 enzyme inhibitory capacity (1A2, 2C9, 2C19, 2D6, 3A4 subtypes), skin permeation (log Kp) that obtained from Swiss ADME prediction are listed in Table 3. Secondly, Discovery Studio software package was used to calculate the ADME/Tox properties. The results of Lipinski’s rule calculation for compounds 1, 2 and 3, include MW, number of HBA, number of HBD, polar surface area (PSA), QP log S that obtained from Discovery Studio prediction are shown in Table 4. The ADME prediction results for compounds 1, 2 and 3, cytochrome P450 2D6 (CYP2D6) enzyme inhibitory capacity, hepatotoxicity, plasma protein binding (PPB) ability that obtained from Discovery Studio prediction are showed in Table 5. The Swiss ADME prediction results showed that the molecular weight values of the compounds were more than 500 (should be between 150 and 500 Daltons); the fraction of carbons in the sp3 hybridization for compound 3 was 0.17 (should be not less than 0.25); the number of rotatable bonds for compounds 2 and 3 were 13 and 14, respectively (should be no more than 9); the number of HBA for the compounds was more than 10, and the number of HBD was more than 5; the TPSA values for the compounds were all more than 130 Å2, the XLOGP3 values for compounds 1 and 2 were -0.88 and -1.60 respectively (should be between -0.7 and +5.0). The gastrointestinal (GI) absorption abilities for the compounds were low; the compounds were unable to permeate BBB; compound 1 was substrate of P-glycoprotein (P-gp), and compounds 2 and 3 were not substrate of P-gp; none of the screened compounds could inhibit the cytochrome P450 enzyme, including 1A2, 2C9, 2C19, 2D6, 3A4

subtypes. The Discovery Studio prediction results were compared with the above prediction, and the two prediction results are basically the same in the following respects: the molecular weight values, the number of HBA and HBD, the PSA values, log S, the cytochrome P450 2D6 enzyme inhibitory capacity. It’s worth noting that the unique results calculated from the Discovery Studio prediction: the hepatotoxicity values of all the three compounds were less than -0.4095, indicating they had less or no hepatotoxicity; and the compounds did not have the plasma protein biding (PPB) ability. The screened compounds came from natural compounds, the results of Lipinski’s calculation and ADME prediction indicated that these compounds need further structural optimization. 3.3. Molecular dynamics simulations analysis Molecular dynamics simulations were performed in 100 ns to analyze the steady nature and conformations stability of the PA-ligand complexes (ligands: compounds 1, 2, and 3). The RMSD was used to estimate the stability of a protein-ligand and protein systems. The RMSD trajectories of PA subunit (free protein) and PA subunit-ligand complexes during 100 ns simulation are shown in Fig. 4. The RMSD values of alpha-carbons for both free and bounded protein complexes changed during the first 50 ns simulation, and became stable after 50 ns. For PA-ligand (compounds 1, 2 and 3), there were some fluctuations in the beginning especially compounds 1 and 3 and then all the complexes gradually reached equilibrium at 75 ns simulation. During the first 35 ns simulation, the fluctuation curve of compound 1 was lower than the others, and the curve of compound 2 became lower than the others after 50 ns simulation. Overall, the RMSD curves for compound 2 was more stable than the others, indicating compound 2 could form more stable complex with PA subunit. The torsional conformations of each rotatable bond in the ligand throughout the simulation trajectories were calculated to understand the conformational evolution of the identified ligands (Fig. 5, compound 1; Fig. 6, compound 2; Fig. 7, compound 3). The dartboard plots on the left illustrated the angle of each bond at a given time during simulations. The beginning of the simulations was placed in the center of the radial plot, and the time evolution was plotted radially outwards. The histograms on the right showed the probability of the torsions as a function of angle, where the radial coordinate was the torsion potential of the rotatable bond, and the angular coordinate was the torsional angle (Zhang et al., 2017).

Table 2 The result of Lipinski’s rule calculation for compounds 1, 2 and 3 obtained from Swiss ADME prediction. No.

MW

Fraction Csp3

No. rotatable bonds

No. HBA

No. HBD

TPSA (Å2)

XLOGP3

Log S

1 2 3

610.52 770.73 718.61

0.44 0.57 0.17

7 13 14

16 19 16

10 11 9

269.43 304.21 278.04

−0.88 −1.60 3.98

−2.88 −2.92 −6.22

MW: molecular weight; HBA: H-bond acceptors; HBD: H-bond donors; TPSA: topological polar surface area; XLOGP3: Log Po/w (XLOGP3); Log S: Log S (ESOL). 5

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al.

Table 3 The ADME prediction results for compounds 1, 2 and 3 obtained from Swiss ADME prediction. No.

Solubility (mg/mL)

GI absorption

BBB permeant

P-gp substrate

CYP1A2 inhibitor

CYP2C19 inhibitor

CYP2C9 inhibitor

CYP2D6 inhibitor

CYP3A4 inhibitor

Log Kp (cm/s)

1 2 3

7.97e-01 9.33e-01 4.33e-04

Low Low Low

No No No

Yes No No

No No No

No No No

No No No

No No No

No No No

−10.65 −12.14 −7.86

GI: gastrointestinal; BBB: blood brain barrier; P-gp: P-glycoprotein; Log Kp: represented skin permeation.

groups and residue Lys643, Gly622, Thr639, Glu623 and Ser709. The torsional angle changes in l and m promoted the π-π interaction between the 4H-chromen-4-one and residue Trp706. Similarly, the torsional angle changes in n corresponded to the π-π interaction between the aromatic ring and residue Phe411. The dial plots and the bar charts showed that the bonds l, m and n had stronger rigidity, and bonds a and b were more flexible. Compounds 2 and 3 bind similarly to the PA subunit as shown in Fig. 2b and Fig. 6, Fig. 3c and Fig. 7. Taking compound 2 in PA subunit binding cavity as the example, the torsional angle changes in a and b promoted the hydrogen-bond interaction between the hydroxyl groups and residue Asn412, the torsional angle changes in g, h, o, p and w corresponded to the hydrogen-bond interaction between the hydroxyl groups and residue Glu623, Glu630, Lys643 and Asn703. The dial plots and the bar charts showed that the bond g had stronger rigidity, and bonds b and h were more flexible. In Fig. 2c and Fig. 7, as for compound 3 in PA subunit binding cavity, the torsional angle changes in t promoted the hydrogen-bond interaction between the oxygen atom and residue Lys643, the torsional angle changes in b and m corresponded to the hydrogen-bond interaction between the hydroxyl groups and residue Cys415 and Glu623. The torsional angle changes in q and s promoted the π-π interaction between the aromatic ring and residue Lys643, the torsional angle changes in l corresponded to the π-π interaction between the aromatic ring and residue Phe411. The dial plots and the bar charts showed that the bonds q, s and t had stronger rigidity, and bond m was more flexible. The rotatable bonds and their rigid/flexible were shown in Fig. 5, Fig. 6, and Fig. 7, compound 3 had fewer rotatable bonds, indicating that the molecule had strong rigidity. The complex was not stable after binding to PA subunit. Residues Glu623 and Lys643 played a crucial role in forming water bridges and hydrogen bonds with compound 1, the compound and Lys643, Asn703 formed water bridges and hydrogen bonds, maintained effectively by 108 % and 116 % respectively of the interactions fraction. In Fig. 8 for compound 1, these stronger interactions were also clearly observed: the compound and Thr639 formed water bridges; the compound and Thr706 formed mostly hydrophobic interactions, a small part of water bridges and hydrogen-bond interactions; the compound and Ser709 formed water bridges and hydrogen-bond interactions. The residues Glu623 and Lys643 of compound 2 formed water bridges and hydrogen bonds, maintained effectively by 80 % and 88 % respectively of the interactions fraction; the compound and Asn703 formed mostly water bridges, maintained effectively by 18 % of the interactions fraction with a weak interaction. In Fig. 8, the compound and Glu630 formed mostly hydrogen bond interactions, a small part of water bridges was also clearly observed. The Glu623 of compound 3 formed water bridges and hydrogen bonds, maintained effectively by 96 % of the interactions fraction; the compound and Lys643 formed mostly water bridges, a small part of hydrogen bond interactions, maintained effectively by 30 % of the interactions fraction with a weak interaction; and the compound and Asn703 formed water bridges, maintained effectively by 50 % of the interactions fraction. In Fig. 8, the compound and Gln408 formed hydrophobic interactions was also clearly observed.

Table 4 The result of Lipinski’s rule calculation for compounds 1, 2 and 3 obtained from Discovery Studio prediction. No.

MW

No. HBA

No. HBD

PSA (Å2)

QP log S

1 2 3

610.524 770.737 716.608

21.5 27.1 22

9 11 7

264.908 301.26 310.894

−1.764 −3.336 −3.724

MW: molecular weight; HBA: H-bond acceptors; HBD: H-bond donors; PSA: polar surface area. Table 5 The ADME prediction results for compounds 1, 2 and 3 obtained from Discovery Studio prediction. No.

CYP2D6 Prediction

Hepatotoxic

PPB Prediction

1 2 3

False False False

−6.60977 −8.01003 −7.14428

False False False

EXT_CYP2D6_Prediction: false: non-inhibitor; true: inhibitor. EXT_ Hepatotoxic: < -0.4095: non-toxic; > -0.4095: inhibitor. EXT_PPB_Prediction: plasma protein binding ability, false: ≥ 90 %; true: ≤ 90 %.

Fig. 4. The RMSD trajectories of PA subunit-ligand complexes during 100 ns simulations.

The hydrogen-bond and other interactions are directional, the interaction between receptor and ligand must satisfy the conditions of distance and direction. The relationships between the torsional angle changes and the interactions were discovered by combining Fig. 2, Fig. 5, Fig. 6 and Fig. 7. In Fig. 2a and Fig. 5, as for compound 1 in PA subunit binding cavity, the torsional angle changes in a and b promoted the hydrogenbond interaction between the hydroxyl groups and residue Asn703. In a similar way, the torsional angle changes in c, j, o, p and q were corresponding to the hydrogen-bond interaction between the hydroxyl 6

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al.

Fig. 5. Graphics of the torsional conformation of each rotatable bond of the ligand compound 1 in PA subunit pocket during MD simulations. The schematic 2D diagram of each ligand is shown with color-coded rotatable bonds on the top. The dial plots of the angle of each bond after simulations are displayed on the left, and the bar charts of the torsional probability as a function of angle are shown on the right.

Based on the above analysis, compound 1 formed more stable interactions than compounds 2 and 3. Comparing Fig. 2 with Fig. 8, the important hydrogen bonds which have been shown in the initial docking models did not change during the MD simulations. The hydrogen-bond interactions between the ligands and residues Glu623, Lys643, Asn703 were also observed in Fig. 8. The water bridge has also been observed as an important interaction between the ligands and above residues. The hydrophobic interactions between the ligands and residues Trp706 were should not be ignored in Fig. 8. Hence, that four residues Glu623, Lys643, Asn703 and Trp706 are identified as the key residues. The result of the calculation is in accordance with the experimental conclusion from the literatures (Massari et al., 2013; Pagano et al., 2013; Tintori et al., 2014; Trist et al., 2016). Six properties were examined to illustrate the stabilities of the selected ligands in the PA subunit binding pocket during the simulation of 100 ns as shown in Fig. 9 (a), (b) and (c): (1) Ligand RMSD: Root mean square deviation of a ligand with respect to the reference conformation (typically the first frame is used as the reference and it is regared as time t = 0); (2) Radius of gyration (rGyr): It is used to measure the ‘extendedness’ of a ligand, and is equivalent to its principal moment of inertia; (3) intramolecular hydrogen bond (intraHB): Number of internal hydrogen bonds (HB) within a ligand molecule. (4) Molecular surface area (MolSA): Molecular surface was calculated with 1.4 Å probe radius; (5) Solvent accessible surface area (SASA): Molecular surface area of accessible by a water molecule; (6) Polar surface area

(PSA): Solvent accessible surface area in a molecule contributed only by oxygen and nitrogen atoms (Raj et al., 2015). As shown in the Fig. 9 (a), (b) and (c), the RMSD values of compound 1 were below 2 Å; the RMSD values of compound 2 were higher than 2 Å, and remained generally stable; the values of compound 3 were nearly 3 Å from 0 ns to 50 ns, and the values were higher than 3 Å after 50 ns. The rGyr values for the ligand compound 1 in the PA subunit binding pocket were 4.4–4.7 Å from 0 ns to 20 ns, and the radius values were 4.2–4.5 Å after 20 ns. The rGyr values for compound 2 in the same binding pocket were basically stable at 6.0 Å from 0 ns to 40 ns, and after 65 ns; the values were fluctuated between 5.6 and 6.2 Å from 40 ns to 65 ns. The rGyr values for compound 3 could be divided into two distinct stages, the values were 4.6–5.0 Å from 0 ns to 40 ns, and the values were basically stable at 5.1–5.25 Å after 40 ns. The constant values showed steady behavior. The intraHB, MolSA, SASA and PSA plots also indicated the consistency of the ligand during the whole simulation process. For properties of MolSA, SASA and PSA for the compound 1, there were upward trend fluctuation occurred from 30 ns to 50 ns, and lower fluctuation occurred from 50 ns to 70 ns. In the rGyr, MolSA, SASA and PSA plots for the compound 3, there was obvious shaky at 20 ns consistently. The RMSD, rGyr and PSA curves were smooth after 50 ns. These curves discovered the consistency of the ligands in the PA subunit binding pocket over the simulation trajectory.

7

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al.

Fig. 6. Graphics of the torsional conformation of each rotatable bond of the ligand compound 2 in PA subunit pocket during MD simulations.

compounds 1-PA subunit complex, -58.18 kcal/mol for compounds 2PA subunit complex, and -68.64 kcal/mol for compounds 3-PA subunit complex. Out of the three compounds, compound 1 showed the lowest total binding free energy. Among all the interactions, the contribution of the covalent binding and the van der Waals energy were more than other items. The values of H-bond interaction for compounds 1, 2, 3 were -4.64, -4.20, -4.21 kcal/mol, respectively, it was obvious that compound 1 formed more stable hydrogen bonds with the residues of PA subunit binding pocket.

3.4. Calculation for prime molecular mechanics / generalized born surface area (MM-GBSA) The total binding free energy was contributed from the Coulomb energy (Coulomb), the covalent binding (Covalent), the hydrogen bonding (Hbond), the lipophillic binding (Lipo), the π-π packing interaction (Packing), the solvent generalized binding and the binding from the van der Waals energy (VDW). The contribution of each item was indicated in Fig. 10. The total binding free energy values between compounds 1, 2, 3 and PA subunit binding pocket were calculated by MM-GBSA method, respectively. The binding free energy value was -77.00 kcal/mol for 8

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al.

Fig. 7. Graphics of the torsional conformation of each rotatable bond of the ligand compound 3 in PA subunit pocket during MD simulations.

4. Conclusion

binding free energy, forming more stable hydrogen bonds with the residues of PA subunit binding pocket. Besides, key amino acid residues of IAV RNA polymerase PA subunit were identified, which would provide a good basis for further research on novel and effective influenza virus PA-PB1 interaction inhibitors. Considering the fact that no RNA polymerase inhibitors in market so far, our research may help to the deepen understanding the inhibitory mechanism of RNA polymerase at atomic level, as well as to expand the chemical space of RNA polymerase inhibitors.

Three natural compounds were identified as potential inhibitors for the influenza virus PA-PB1 subunits interaction from the natural resource. The stability of the complexes of PA-hits was validated through molecular dynamics (MD) simulation, and the protein–ligand contacts over the course of trajectory showed four key amino acid residues (Glu623, Lys643, Asn703 and Trp706) of IAV. The binding free energy and contributed energies were calculated by using MM-GBSA method. Out of the three compounds, compound 1 showed the lowest total

9

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al.

Fig. 8. The bar charts of protein–ligand contacts and contacts over the course of trajectory. The legend for each color section is given.

Author statement

There are no conflicts of interests and the contents of this manuscript has not been published or submitted to any journal aside from your journal. We will be happy to receive any questions or comment from you and please don't hesitate to contact us anytime for clarification.

We would like to re-submit the revised manuscript entitled “Structure-based virtual screening of influenza virus RNA polymerase inhibitors from natural compounds: molecular dynamics simulation and MM-GBSA calculation”, for publication in Computational Biology and Chemistry as original article. On checking our manuscript CBAC_2019_864_R1 entitled "Structurebased virtual screening of influenza virus RNA polymerase inhibitors from natural compounds: molecular dynamics simulation and MMGBSA calculation", we have revised the manuscript according to the referee’s comments.

Declaration of Competing Interest All authors declare that there are no conflicts of interest associated with the publication of this manuscript.

10

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al.

Fig. 10. MM-GBSA method was used to predict the binding free energy between the PA subunit and the ligands (compounds 1, 2 and 3). The total MMGBSA binding free energy was contributed from the coulomb energy (Coulomb), the covalent binding (Covalent), the hydrogen bonding (Hbond), the lipophillic binding (Lipo), the π-π packing interaction (Packing), the solvent generalized binding and the binding from the van der Waals energy (VDW).

Acknowledgements The work was supported by the Fund for the Natural Science Foundation of Liaoning province (Grant No. 20170540854), the Natural Science Foundation of Liaoning province (Grant No. 201602695 and 2019-ZD-0458), the Scientific Research Foundation of Educational Department of Liaoning Province (Grant No. L2015517), the Innovation and Entrepreneurship Training Program for College Students of Shenyang Pharmaceutical University (Grant No. 201910163044). Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.compbiolchem.2020. 107241. References Aier, I., Varadwaj, P., Raj, U., 2016. Structural insights into conformational stability of both wild-type and mutant EZH2 receptor. Sci. Rep. 6, 34984. Bouvier, N., Palese, P., 2008. The biology of influenza viruses. Vaccine 26 (Suppl 4), D49–D53. Chang, S., Park, J., Kim, Y., Kang, J., Min, J.-Y., 2016. A natural component from Euphorbia humifusa Willd displays novel, broad-spectrum anti-influenza activity by blocking nuclear export of viral ribonucleoprotein. Biochem. Biophys. Res. Commun. 471 (2), 282–289. Chen, J., Jiang, H., Li, F., Hu, B., Wang, Y., Wang, M., Wang, J., Cheng, M., 2018. Computational insight into dengue virus NS2B-NS3 protease inhibition: a combined ligand- and structure-based approach. Comput. Biol. Chem. 77, 261–271. Daina, A., Michielin, O., Zoete, V., 2017. SwissADME: a free web tool to evaluate pharmacokinetics, druglikeness and medicinal chemistry friendliness of small molecules. Sci. Rep. 7, 42717. Deyde, V., Xu, X., Bright, R., Shaw, M., Smith, C., Zhang, Y., Shu, Y., Gubareva, L., Cox, N., Klimov, A., 2007. Surveillance of resistance to adamantanes among influenza A(H3N2) and A(H1N1) viruses isolated worldwide. J. Infect. Dis. 196 (2), 249–257. Dias, A., Bouvier, D., Crépin, T., McCarthy, A., Hart, D., Baudin, F., Cusack, S., Ruigrok, R., 2009. The cap-snatching endonuclease of influenza virus polymerase resides in the PA subunit. Nature 458 (7240), 914–918. Friesner, R., Murphy, R., Repasky, M., Frye, L., Greenwood, J., Halgren, T., Sanschagrin, P., Mainz, D., 2006. Extra precision Glide: docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes. J. Med. Chem. 49 (21), 6177–6196. He, X., Zhou, J., Bartlam, M., Zhang, R., Ma, J., Lou, Z., Li, X., Li, J., Joachimiak, A., Zeng, Z., Ge, R., Rao, Z., Liu, Y., 2008. Crystal structure of the polymerase PAC-PB1N complex from an avian influenza H5N1 virus. Nature 454, 1123–1126. Hensens, O., Goetz, M., Liesch, J., Zink, D., Raghoobar, S., Helms, G., Singh, S., 1995. Isolation and structure of flutimide, a novel endonuclease inhibitor of influenza virus. Tetrahedron Lett. 36 (12), 2005–2008. Hilleman, M., 2002. Realities and enigmas of human viral influenza: pathogenesis,

Fig. 9. Variation in the ligands (compounds 1, 2 and 3) properties w.r.t time during the course of 100 ns simulation.

11

Computational Biology and Chemistry 85 (2020) 107241

Z. Jin, et al. epidemiology and control. Vaccine 20, 3068–3087. Itteboina, R., Ballu, S., Sivan, S., Manga, V., 2017. Molecular modeling-driven approach for identification of Janus kinase 1 inhibitors through 3D-QSAR, docking and molecular dynamics simulations. J. Recept. Signal Transduct. 37 (5), 453–469. Kaplan, M., 2010. The 2009 H1N1 swine Flu pandemic: reconciling goals of patents and public health initiatives. 20 Fordham Intell. Prop. Media & Ent. L.J. 991. Li, R., Liu, T., Liu, M., Chen, F., Liu, S., Yang, J., 2017b. Anti-influenza A virus activity of dendrobine and its mechanism of action. J. Agric. Food Chem. 65 (18), 3665–3674. Liu, M., Lam, M., Zhang, Q., Elderfield, R., Barclay, W., Shaw, P.-C., 2015. The functional study of the N-terminal region of influenza B virus nucleoprotein. PLoS One 10 (9), e0137802/1–e0137802/12. Liu, M., Lo, C.-Y., Wang, G., Chow, H.-F., Ngo, J., Wan, D., Poon, L., Shaw, P.-C., 2017a. Identification of influenza polymerase inhibitors targeting polymerase PB2 capbinding domain through virtual screening. Antiviral Res. 144, 186–195. Liu, Y.-Y., Feng, X.-Y., Jia, W.-Q., Jing, Z., Xu, W.-R., Cheng, X.-C., 2019. Identification of novel PI3Kδ inhibitors by docking, ADMET prediction and molecular dynamics simulations. Comput. Biol. Chem. 78, 190–204. Massari, S., Nannetti, G., Goracci, L., Sancineto, L., Muratore, G., Sabatini, S., Manfroni, G., Mercorelli, B., Cecchetti, V., Facchini, M., Palù, G., Cruciani, G., Loregian, A., Tabarrini, O., 2013. Structural investigation of cycloheptathiophene-3-carboxamide derivatives targeting influenza virus polymerase assembly. J. Med. Chem. 56 (24), 10118–10131. Ng, A., Wang, J.-H., Shaw, P.-C., 2009. Structure and sequence analysis of influenza A virus nucleoprotein. Sci. China Ser. C-life Sci. 52 (5), 439–449. Obayashi, E., Yoshida, H., Kawai, F., Shibayama, N., Kawaguchi, A., Nagata, K., Tame, J., Park, S.-Y., 2008. The structural basis for an essential subunit interaction in influenza virus RNA polymerase. Nature 454, 1127–1132. Pagano, M., Castagnolo, D., Bernardini, M., Fallacara, A.L., Laurenzana, I., Deodato, D., Kessler, U., Beatrice Pilger, B., Stergiou, L., Strunze, S., Tintori, C., Botta, M., 2013. The fight against the influenza A virus H1N1: synthesis, molecular modeling, and biological evaluation of benzofurazan derivatives as viral RNA polymerase inhibitors. ChemMedChem 9 (1), 129–150. Raj, U., Kumar, H., Gupta, S., Varadwaj, P., 2015. Novel DOT1L receptornatural inhibitors involved in mixed lineage leukemia: a virtual screening, molecular docking and dynamics simulation study. Asian Pacific J. Cancer Prev. 16 (9), 3817–3825. Shivakumar, D., Williams, J., Wu, Y., Damm, W., Shelley, J., Sherman, W., 2010. Prediction of absolute solvation free energies using molecular dynamics free energy perturbation and the OPLS force field. J. Chem. Theory Comput. 6 (5), 1509–1519. Singh, S., Tomassini, J., 2001. Synthesis of natural flutimide and analogous fully substituted pyrazine-2,6-diones, endonuclease inhibitors of influenza virus. J. Org. Chem. 66, 5504–5516. Sugita, Y., Noda, T., Sagara, H., Kawaoka, Y., 2011. Ultracentrifugation deforms unfixed influenza A virions. J. Gen. Virol. 92, 2485–2493. Suzuki, Y., 2005. Sialobiology of influenza: molecular mechanism of host range variation

of influenza viruses. Biol. Pharm. Bull. 28 (3), 399–408. Tintori, C., Laurenzana, I., Fallacara, A., Kessler, U., Pilger, B., Stergiou, L., Botta, M., 2014. High-throughput docking for the identification of new influenza A virus polymerase inhibitors targeting the PA–PB1 protein–protein interaction. Bioorg. Med. Chem. Lett. 24 (1), 280–282. Tong, S., Zhu, X., Li, Y., Shi, M., Zhang, J., Bourgeois, M., Yang, H., Chen, X., Recuenco, S., Gomez, J., Chen, L.-M., Johnson, A., Tao, Y., Dreyfus, C., Yu, W., McBride, R., Carney, P., Gilbert, A., Chang, J., Guo, Z., Davis, C., Paulson, J., Stevens, J., Rupprecht, C., Holmes, E., Wilson, I., Donis, R., 2013. New world bats harbor diverse influenza A viruses. PLoS Pathog. 9 (10), e1003657. Trist, I., Nannetti, G., Tintori, C., Fallacara, A., Deodato, D., Mercorelli, B., Palù, G., Wijtmans, M., Gospodova, T., Edink, E., Verheij, M., De Esch, I., Viteva, L., Loregian, A., Botta, M., 2016. 4,6-Diphenylpyridines as promising novel anti-influenza agents targeting the PA-PB1 protein–protein interaction: structure-activity relationships exploration with the aid of molecular modeling. J. Med. Chem. 59 (6), 2688–2703. van de Sandt, C., Bodewes, R., Rimmelzwaan, G., de Vries, R., 2015. Influenza B viruses: not to be discounted. Future Microbiol. 10 (9), 1447–1465. Wang, F.-F., Yang, W., Shi, Y.-H., Le, G.-W., 2016. In silico study on β-aminoketone derivatives as thyroid hormone receptor inhibitors: a combined 3D-QSAR and molecular docking study. J. Biomol. Struct. Dyn. 34 (12), 2619–2631. Wang, M., Li, W., Wang, Y., Song, Y., Wang, J., Cheng, M., 2018a. In silico insight into voltage-gated sodium channel 1.7 inhibition for anti-pain drug discovery. J. Mol. Graph. Model. 84, 18–28. Wang, M., Wang, Y., Kong, D., Jiang, H., Wang, J., Cheng, M., 2018b. In silico exploration of aryl sulfonamide analogs as voltage-gated sodium channel 1.7 inhibitors by using 3D-QSAR, molecular docking study, and molecular dynamics simulations. Comput. Biol. Chem. 77, 214–225. Wang, Y., Feng, S., Gao, H., Wang, J., 2019. Computational investigations of gram-negative bacteria phosphopantetheine adenylyltransferase inhibitors using 3D-QSAR, molecular docking and molecular dynamic simulations. J. Biomol. Struct. Dyn. https://doi.org/10.1080/07391102.2019.1608305. Wu, Y., Wu, Y., Tefsen, B., Shi, Y., Gao, G., 2014. Bat-derived influenza-like viruses H17N10 and H18N11. Trends Microbiol. 22 (4), 183–191. Yang, C.-H., Tan, D.-H., Hsu, W.-L., Jong, T.-T., Wen, C.-L., Hsu, S.-L., Chang, P.-C., 2014. Anti-influenza virus activity of the ethanolic extract from Peperomia sui. J. Ethnopharmacol. 155 (1), 320–325. Zhang, J., Liu, X., Wang, S.-Q., Liu, G.-Y., Xu, W.-R., Cheng, X.-C., Wang, R.-L., 2017. Identification of dual ligands targeting angiotensin II type 1 receptor and peroxisome proliferator-activated receptor-γ by core hopping of telmisartan. J. Biomol. Struct. Dyn. 35 (12), 2665–2680. Zhou, J., Zou, L., Zhang, X., Liao, J., Ni, H., Hou, N., Wang, Y., Li, H., Wu, J., Jonges, M., Meijer, A., Koopmans, M., Ke, C., 2011. Adamantone- and Oseltomivir- resistant seasonal A (H1N1) and pandemic A (H1N1) 2009 influenza viruses in Guangdong, China, during 2008 and 2009. J. Clin. Microbiol. 49 (7), 2651–2655.

12