Molecular Docking and Virtual Screening of an Influenza Virus Inhibitor That Disrupts Protein–Protein Interactions

Influenza is an acute respiratory infection caused by the influenza virus, but few drugs are available for its treatment. Consequently, researchers have been engaged in efforts to discover new antiviral mechanisms that can lay the foundation for novel anti-influenza drugs. The viral RNA-dependent RNA polymerase (RdRp) is an enzyme that plays an indispensable role in the viral infection process, which is directly linked to the survival of the virus. Methods of inhibiting PB1–PB2 (basic polymerase 1–basic polymerase 2) interactions, which are a key part of RdRp enzyme activity, are integral in the design of novel antiviral drugs, a specific PB1–PB2 interactions inhibitor has not been reported. We have screened Enamine’s database and conducted a parallel screening of multiple docking schemes, followed by simulations of molecular dynamics to determine the structure of a stable ligand—PB1 complex. We also calculated the free energy of binding between the screened compounds and PB1 protein. Ultimately, we screened and identified a potential PB1–PB2 inhibitor using the ADMET prediction model.


Introduction
The influenza virus causes acute respiratory viral infections, which induce complications that often require hospitalization and may even lead to death. At the same time, the virus is easily transmissible from person to person and can infect individuals of all ages [1]. Consequently, seasonal influenza is deemed a serious public health concern.
Influenza viruses belong to the family of Orthomyxoviridae and are classified into four types: A, B, C, and D [2][3][4], and the first three types are able to infect humans. Influenza A virus (IAV) is prone to cause periodic pandemics because of the frequent mutation to escape the host immune system [5,6]. Influenza B virus (IBV) is relatively pathogenic to human beings, but it has not caused a global pandemic [2]. Influenza C virus (ICV) only causes mild infection, and affects public health slightly [7]. Influenza D virus (IDV) mainly infects pigs and cattle, but not human [8].
Current treatment options for tackling influenza viruses are limited, and drug resistance is a growing problem. The WHO recommends treatment with a neuraminidase inhibitor (oseltamivir) within 48 h of flu symptoms [9]. All circulating influenza viruses are resistant to adamantane drugs, such as amantadine and rimantadine [2,10,11]. According to the findings of the Chinese National Influenza Center's resistance surveillance, all strains of influenza A/H1N1 and A/H3N2 viruses subtypes are resistant to amantadine analogues [10,11], and resistance to neuraminidase inhibitors has been reported in three cases, but the surveillance data revealed one case of resistance to neuraminidase inhibitors among influenza A/H1N1 viruses strains between 1 October 2019 and 13 September 2020 [12]. Therefore, the development of novel anti-influenza drugs is an urgent task.
Up to now, the antiviral drugs in the international market have been mainly divided into three categories: M2 ion channel blockers, neuraminidase inhibitors and RdRp drugs. Amantadine and Rimantadine are M2 ion channel blockers which block ion transport by binding to M2 ion channel proteins, thereby inhibiting viral replication and preventing the viral infection of new cells. Due to the S31N mutation [13,14], Amantadine cannot tightly fit with the ion channel, allowing the ions to pass through, and finally leading to drug resistance.
The neuraminidase inhibitors mainly include Zanamivir, Oseltamivir and Peramivir, which can mimic the natural substrate of NA, sialic acid, by blocking its active site and preventing it from catalyzing the hydrolysis of sialic acid, thus preventing the release of viral particles. Drug resistance is mainly due to the mutation of the binding site and the reduction of neuraminidase activity [15,16].
Recently developed anti-influenza drugs are mainly RdRp drugs that target the virus [17,18]. RdRp of influenza A virus consists of three different subunits: acidic polymerase, basic polymerase 1, and basic polymerase 2 (PA, BP1, and PB2), and antiviral effects can be achieved through the successful blocking of protein-protein interactions (PPIs) during protein assembly [19,20]. The currently marketed RdRp drugs include Favipiravir and Baloxavir. Favipiravir, as nucleoside analogues, is only used after the conventional antiviral drug treatment is ineffective, owing to its disadvantages of high toxicity, teratogenicity and abnormal behavior, etc. Baloxavir is an endonuclease inhibitor superior to neuraminidase inhibitors with a lower rate of adverse reactions. Reports indicate that the I38T/M/F mutation may lead to the development of resistance to Baloxavir. In recent years, small molecule inhibitors targeting PA-PB1 interactions have been extensively studied [21][22][23]. However, inhibitors targeting PB1-PB2 interactions have been rarely reported [18].
In this study, we developed a novel PPI inhibitor via a computer-aided design based on PB1-PP2 interactions [24,25]. Figure 1 depicts the positions of interacting PB1 and PB2 (the red part is PB1, and the cyan part is PB2). The end portion of PB2 is wrapped in three α-helices at the C-terminal of PB1. An antiviral effect can be expected through the occupation of the PB1 pocket by small molecules that prevent the binding of PB1 and PB2. To disrupt the interactions between PB1 and PB2, we examined key amino acids at their junction. For the PB1 fraction, Pooleal et al. characterized a polypeptide from the Cterminus of PB1, which disrupted the binding of PB1C to PB2N in experiments, thus inhibiting the activity of influenza RdRp and effectively inhibiting virus replication [26]. Monoclonal antibodies against the PB1 polypeptide of the influenza A virus can also interfere with viral mRNA synthesis in vitro. Moreover, the deletion and mutation of PB1 binding site residues can disrupt normal viral replication [27,28]. To disrupt the interactions between PB1 and PB2, we examined key amino acids at their junction. For the PB1 fraction, Pooleal et al. characterized a polypeptide from the C-terminus of PB1, which disrupted the binding of PB1C to PB2N in experiments, thus inhibiting the activity of influenza RdRp and effectively inhibiting virus replication [26]. Monoclonal antibodies against the PB1 polypeptide of the influenza A virus can also interfere with viral mRNA synthesis in vitro. Moreover, the deletion and mutation of PB1 binding site residues can disrupt normal viral replication [27,28].
For the PB2 fraction, viral RNA polymerase activity was decreased through the deletion of 27 amino acids from the PB2N terminal [26,29]. These 27 amino acids, which are located at the PB2N terminal, can be applied as a chromogen reagent of viral transcription Viruses 2021, 13, 2229 3 of 14 after binding to green fluorescent protein [30]. A monoclonal antibody developed against the influenza A virus PB2 peptide also interferes with PB1-PP2 binding, thereby inhibiting viral activity in vitro [31,32].
There is a high level of recognition specificity related to PB1 and PB2, and viral replication is inhibited by mutual loss, mutation, and substitution. Therefore, the disruption of the binding of these parts (PB1 731-757 and PB2 0-12) can effectively inhibit viral replication. We therefore designed small molecule inhibitors by replacing PB2.
In this study, we have designed a novel small molecule PPI inhibitor to occupy the PB1-PB2 interactions pocket with the aid of a computer design program. We subsequently conducted the virtual screening of Enamine PPI library containing 40,640 compounds, ultimately identified a molecule that could potentially serve as a lead compound.

Pretreatment of the PB1 and PB2 Proteins of the Influenza Virus
The X-ray crystallographic structure of RNA polymerase PB1-PB2 subunits [34] (PDB ID: 2ZTT, resolution: 2.10 Å, r-value free: 0.272) was downloaded from the Protein Data Bank [35]. We retained one of the PB1 chains for docking and replaced selenomethionine with methionine in the PDB file format. The further processing of proteins was conducted according to docking software requirements. Using the AutoDock Vina software (Scripps Research, San Diego, CA, USA) [36], we manually removed water molecules, ligands, and ions, hydrogenated and protonated each residue in the protein file using PyMOL [37], finally converted the PDB to a PDBPT file for subsequent docking using the Obabel program (University of Pittsburgh, Pittsburgh, USA) [38]. We applied the existing pre-processing function in the SYBYL software to remove water molecules, ligands, and ions from the protein file and to hydrogenate each residue. We also used this to repair end residues, reinsert position files, repair residues, and minimize energy and possible interatomic collisions. We used the protein preparation tool in the Discovery Studio software to complete the process of removing water molecules, ligands, and ions, hydrogenated and protonated each residue in the protein file.

Ligand Preparation
We selected Enamine's database, which provides a specific molecular library for PPI screening [39]. This database was selected because the molecules within this collection have an α-helix-like structure, with the ability to mimic protein motifs. The compounds in the database can be purchased directly for testing after completion of the screening process. The molecular library was split using the OpenBabel program and was subsequently hydrogenated, and energy minimized using a high-precision MMFF94 force field [38]. The 3D structures ultimately generated were used in the molecular docking studies.

Molecular Docking
To select the docking pocket ( Figure 2), part of the PB2 short peptide (0-12), which plays a major role in PPI, was selected as the simulated binding site using the SYBYL and Discovery Studio programs, and the amino acid space within 5 Å of the short peptide was selected as the docking pocket. A cube centered at x = 4.639, y = 7.639, and z = 5.861 with the side lengths of x, y, and z set respectively at 18.  The database molecules were docked into the docking pockets of PB1 protein using the Surflex-Dock Screen provided in the Sybyl 2.0 software. The docking results were scored and ranked, and 2000 compound conformations with high scores were selected for further docking. These molecular conformations were docked again to the active site by applying three modes of docking: the Surflex-Dock GeomX mode in the SYBYL 2.0 software, the LibDock mode in the Discovery Studio software (3dexperiencem, Aachen, Germany), and the AutoDock Vina software. Three molecules were selected for further study according to their docking positions relating to the receptor, considered in combination with their scores.

Simulation of Molecular Dynamics
Molecular dynamics (MD) simulations were performed using the docking results to enable the further analysis of the binding of protein-ligand complexes. All MD simulations were performed using the GROMACS package (University of Groningen, Groningen, Netherlands), version 2020.3 [40][41][42]. The AMBER99SB-ILDN force field was used to optimize proteins [43], and water molecules were modeled with TIP3P [44,45]. We used ACPYPE to calculate ligand charges and generate topology files for GAFF force [46]. The MD box was set up as a cube with the distance between each atom and the box exceeding 0.8 nm. SOL water was added to the box at a density of 1000 g/L. Chloride ions were used to replace water molecules in the system to achieve an electrically neutral simulation system. The relaxation system was then simulated using the Verlet algorithm for energy minimization to eliminate the overlap between atoms and avoid problems such as the overly close proximity of homogeneous charges [47]. The system was subjected to a 100 ps restricted kinetic simulation at 298.15 K. Lastly, the regular dynamics of the system were simulated using the Verlet algorithm [47]. The integration step was set at 0.002 ps, and simulations were conducted for 10 ns in total. The simulations were performed in an isothermal isobaric regime at 298.15 K and under 1 bar pressure, with temperature and pressure respectively controlled with the V-rescale and Parrinello-Rahman methods [48], and PBC (periodic boundary condition) was enabled. The RMSD (root mean squared deviation) was calculated for protein-protein and protein-small molecules interactions. MD trajectories were viewed using VMD software (University of Illinois at Urbana-Champaign, Urbana, IL, USA) [49]. The database molecules were docked into the docking pockets of PB1 protein using the Surflex-Dock Screen provided in the Sybyl 2.0 software. The docking results were scored and ranked, and 2000 compound conformations with high scores were selected for further docking. These molecular conformations were docked again to the active site by applying three modes of docking: the Surflex-Dock GeomX mode in the SYBYL 2.0 software, the LibDock mode in the Discovery Studio software (3dexperiencem, Aachen, Germany), and the AutoDock Vina software. Three molecules were selected for further study according to their docking positions relating to the receptor, considered in combination with their scores.

Simulation of Molecular Dynamics
Molecular dynamics (MD) simulations were performed using the docking results to enable the further analysis of the binding of protein-ligand complexes. All MD simulations were performed using the GROMACS package (University of Groningen, Groningen, Netherlands), version 2020.3 [40][41][42]. The AMBER99SB-ILDN force field was used to optimize proteins [43], and water molecules were modeled with TIP3P [44,45]. We used ACPYPE to calculate ligand charges and generate topology files for GAFF force [46]. The MD box was set up as a cube with the distance between each atom and the box exceeding 0.8 nm. SOL water was added to the box at a density of 1000 g/L. Chloride ions were used to replace water molecules in the system to achieve an electrically neutral simulation system. The relaxation system was then simulated using the Verlet algorithm for energy minimization to eliminate the overlap between atoms and avoid problems such as the overly close proximity of homogeneous charges [47]. The system was subjected to a 100 ps restricted kinetic simulation at 298.15 K. Lastly, the regular dynamics of the system were simulated using the Verlet algorithm [47]. The integration step was set at 0.002 ps, and simulations were conducted for 10 ns in total. The simulations were performed in an isothermal isobaric regime at 298.15 K and under 1 bar pressure, with temperature and pressure respectively controlled with the V-rescale and Parrinello-Rahman methods [48], and PBC (periodic boundary condition) was enabled. The RMSD (root mean squared deviation) was calculated for protein-protein and protein-small molecules interactions. MD trajectories were viewed using VMD software (University of Illinois at Urbana-Champaign, Urbana, IL, USA) [49].

Calculation of the Binding Free Energy
The binding free energy between the protein-ligand complexes was estimated using the MM/PBSA equation [50]. The APBS lattice parameters are output according to the MD results, and the APBS software (Pacific Northwest National Laboratory, Richland, DC, USA) [51] was applied to calculate polar solvation energy (PB) and polar solvation energy (SA). Subsequently, the MM contribution was calculated, and the energy decomposition of the residues was performed. The MM energy was calculated on the basis of directly accumulated potential, taking into account the dielectric constant of the solute (protein). PB energy was calculated by subtracting the energy of the aqueous phase from the energy in the vacuum. SA energy was calculated by multiplying the solvent-accessible surface area of each atom by the surface tension.

Binding Force Analysis
We used the Discovery Studio software to analyze the final results and to examine the non-bonding interactions and hydrogen bonding between the ligand and the receptor.
The non-bonding interactions are the main outlets to examine the hydrophobic interactions between π ring and alkyl. With respect to the geometry of π-π stacked interactions, the distance between the centroid of each π ring pair is less than 6 Å. At least one pair of atoms on both rings is less than 4.5 Å. The angle between the normal of one or both rings and the centroid-centroid vector must be between 0 • and ±50 • , and the angle between the normal to each ring must be between 0 • and ±35 • .
With respect to the geometry of alkyl interactions, the ligand hydrophobic groups are not attached to charged atoms and have a surface area greater than 0.65 times the surface area of the methyl group. Moreover, the distance between the centroid of different hydrophobic groups on the receptor and ligand is less than 4.5 Å.
With respect to the geometry of π-alkyl interactions, the hydrogens acting as the donor must be connected to a nonaromatic carbon atom. The distances between the hydrogen and the center of the π ring must be less than 4.1 Å. The C-H-centroid angle cannot deviate from linearity by more than 20 • and the angle between the C-centroid vector and the normal to the ring plane is not more than 45 • .
With respect to the geometry of hydrogen bond definition, the distance between the hydrogen atom and the hydrogen bond receptor was less than 2.7 or the distance between that hydrogen bond donor and the hydrogen bond receptor is less than 3.4. At the same time, the angle formed between the hydrogen atom-hydrogen bond donor-hydrogen bond receptor should be less than 35 • .

ADMET Prediction
We used the ADMET Descriptors module and the TOPKAT module in Discovery Studio software to predict the ADMET and toxicological properties of the compounds.
The ADMET Descriptors module is developed based on descriptors computed by linear formulas. The water-soluble model was developed using a training set of 775 compounds and a test set of 34 compounds using descriptors of solubility and solubility levels (R 2 = 0.88, SD = 0.79). Blood-brain permeability was developed using a training set of 102 compounds and a test set of 86 compounds. Descriptors used were 2D polar surface area, halogen content, and LogP (R 2 = 0.88, SD = 0.31). A toxicity model was generated by modified Bayesian learning. A leave-one-out cross-validation was conducted. Model sensitivity is 0.816, specificity is 0.854, and concordance rate is 0.834.
Based on the 2D molecular structure, TOPKAT module resolved the predicted compound into functional groups and substructures, and converted them into descriptors, such as two-dimensional molecular, electronic and spatial descriptors, by similarity scoring function. Finally, the descriptors were brought into the different structural toxicity relationship equations of Quan ta va in TOPKAT module to calculate the possible toxicity values of the corresponding compounds.

Protein Preparation
Only 2ZTT and 3A1G have been reported for the crystal structure of RNA polymerase PB1-PB2 subunit [34]. The two protein structures have resolutions between 1.5 Å and 2.5 Å. A Ramachandran plot prediction of the two proteins was generated ( Table 1). Most of the 2ZTT residues were favorably located (96%); a few residues were located in allowed Viruses 2021, 13, 2229 6 of 14 areas (4%), and no residues were outliers. In the case of 3AIG, most of the residues were favorably located (98%), a few residues were in allowed locations (2%), and one residue was an outlier. For this screening, 2ZTT was selected as the docked protein structure. To ensure the generalizability of the screening results, we used the ENDscript/ESPript software to calculate the conservativeness of docked proteins among homologous proteins [52]. As shown in Figure 3, among the homologous proteins, the PB1C-terminal protein of the three α-helices was well conserved sequence comparison has been completed by Clustal Omega web server [53]. The inhibitor designed based on the protein binding site can effectively inhibit different subtypes of influenza A virus and influenza B virus (4WRT, 6F5O and 6QWL). There are some differences (6KUJ) in the corresponding residues of influenza D virus, which may lead to the decrease of activity. The binding site of influenza C virus is too different, so there is no universality.

Protein Preparation
Only 2ZTT and 3A1G have been reported for the crystal structure of RNA polymerase PB1-PB2 subunit [34]. The two protein structures have resolutions between 1.5 Å and 2.5 Å. A Ramachandran plot prediction of the two proteins was generated ( Table 1). Most of the 2ZTT residues were favorably located (96%); a few residues were located in allowed areas (4%), and no residues were outliers. In the case of 3AIG, most of the residues were favorably located (98%), a few residues were in allowed locations (2%), and one residue was an outlier. For this screening, 2ZTT was selected as the docked protein structure. To ensure the generalizability of the screening results, we used the ENDscript/ES-Pript software to calculate the conservativeness of docked proteins among homologous proteins [52]. As shown in Figure 3, among the homologous proteins, the PB1C-terminal protein of the three α-helices was well conserved sequence comparison has been completed by Clustal Omega web server [53]. The inhibitor designed based on the protein binding site can effectively inhibit different subtypes of influenza A virus and influenza B virus (4WRT, 6F5O and 6QWL). There are some differences (6KUJ) in the corresponding residues of influenza D virus, which may lead to the decrease of activity. The binding site of influenza C virus is too different, so there is no universality. Protein pretreatment: given the lack of non-standard amino acid force fields in the docking software, selenomethionine was manually replaced with methionine, and irrelevant atoms were removed during MD simulations and molecular docking. Terminal residues in the protein structure were repaired, hydrogens were added, force field files were reincorporated, and energy minimization of the protein was conducted.

Molecular Docking
Initially, we selected the top 2000 compounds for fine docking, on the basis of the Sybyl docking results. These molecules were docked to the active site using three modes of docking to enable a precise evaluation of the binding. The docking results for the top 200 compounds are shown in Table 2. Protein pretreatment: given the lack of non-standard amino acid force fields in the docking software, selenomethionine was manually replaced with methionine, and irrelevant atoms were removed during MD simulations and molecular docking. Terminal residues in the protein structure were repaired, hydrogens were added, force field files were reincorporated, and energy minimization of the protein was conducted.

Molecular Docking
Initially, we selected the top 2000 compounds for fine docking, on the basis of the Sybyl docking results. These molecules were docked to the active site using three modes of docking to enable a precise evaluation of the binding. The docking results for the top 200 compounds are shown in Table 2. Statistical advantages in docking scores were identified for the following compounds:

Molecular Dynamics
We further checked the conformational plausibility of the screened compounds by performing MD simulations. Specifically, we performed MD simulations for three molecules for 10 ns each. During the simulations, the RMSD values of the entire system showed

Molecular Dynamics
We further checked the conformational plausibility of the screened compounds by performing MD simulations. Specifically, we performed MD simulations for three molecules for 10 ns each. During the simulations, the RMSD values of the entire system showed a gradual convergence, and eventually became stable, as shown in Figure 5.

Molecular Dynamics
We further checked the conformational plausibility of the screened compounds by performing MD simulations. Specifically, we performed MD simulations for three molecules for 10 ns each. During the simulations, the RMSD values of the entire system showed a gradual convergence, and eventually became stable, as shown in Figure 5.

Calculation of the Binding Free Energy
After we obtained stable protein-ligand complexes through MD simulations, we calculated the binding energy between the three compounds and the protein using the MM/PBSA equation (Table 3).

Calculation of the Binding Free Energy
After we obtained stable protein-ligand complexes through MD simulations, we calculated the binding energy between the three compounds and the protein using the MM/PBSA equation (Table 3). As shown in Table 3, the values of the binding free energies of all three systems were below zero, indicating that the binding of protein-ligand complexes was stable. The binding free energies of Compounds 3206 and 4808 were, respectively, −128.049 kJ/mol and −122.761 kJ/mol, which exceeded that of Compound 15144 (−106.93 kJ/mol). Van der Waals interaction accounted for the binding of Compounds 3206 and 4808 with protein, and the electrostatic interaction of these two compounds, had particularly different of −39.663 kJ/mol and 1.721 kJ/mol. These results suggest that compound 4808 is more suitable for use as a lead inhibitor by changing the molecular polarity or local polarity to facilitate increased compound-protein interactions.
The molecular dynamics simulation of the polypeptide-protein complex was performed for 10 ns with the polypeptide from the N-terminus of PB2 as ligands. Considering the macromolecules as the comparison object, the Gibbs free energy was used as the comparison standard. The Gibbs free energy of polypeptide and PB1 was −34.792 kJ/mol comparing with −45.007 kJ/mol of compound 4808 and PB1, which indicated that the screened molecules could inhibit PB1-PB2 interactions effectively and competitively.
Subsequently, we performed a binding free energy decomposition analysis of Compound 4808; the results are depicted in Figure 6.

Binding Force Analysis
In light of the results of the MD simulations, we performed a binding force analysis The results of the analysis of residue decomposition indicated that protein residues 39, 41, and 45 (residue contributions of 0.425, 6.956, and 8.779 kJ/mol, respectively) had an unfavorable impact on the binding. Residues 15, 59, 62, and 66 (residue contributions of −10.113, −8.797, −11.035, and −7.599 kJ/mol, respectively) were the main contributors to the binding of the protein to Compound 4808.

Binding Force Analysis
In light of the results of the MD simulations, we performed a binding force analysis using Discovery Studio's receptor-ligand interaction calculation tool (Figure 7). The interactions between protein and Compound 4808 mainly comprised hydrogen bonding and hydrophobic interactions. For the hydrogen bond interactions, residue 14 formed two carbon-hydrogen bonds with bond lengths of 2.509 Å and 2.534 Å. The following hydrophobic interactions occurred: a π-π stacking interaction between residue 15 and the compound; an alkyl interaction between residue 66 and the compound; and π-alkyl interactions between residues 38, 42, 59, and 62 and the compound.

ADMET Predictions
ADMET predictions indicated that Compound 4808 is slightly soluble in water at 25 °C (log (SW) = −5.137). A brain-blood ratio between 1:1 and 1:5 does not inhibit cytochrome P450 2D6. The Bayesian score for pigment toxicity was −2.63393, which is associated with some toxicity to the liver. The Bayesian score for hepatotoxicity was −1.85674, which is associated with adequate intestinal absorption (ADMET_Absorption_T2_2D <

ADMET Predictions
ADMET predictions indicated that Compound 4808 is slightly soluble in water at 25 • C (log (SW) = −5.137). A brain-blood ratio between 1:1 and 1:5 does not inhibit cytochrome P450 2D6. The Bayesian score for pigment toxicity was −2.63393, which is associated with some toxicity to the liver. The Bayesian score for hepatotoxicity was −1.85674, which is associated with adequate intestinal absorption (ADMET_Absorption_T2_2D < 6.1261).
Toxicity predictions were calculated using the Ames heteroaromatic model. According to the results of the prediction study conducted on PB1-PB2 protein binding, the computed probability of mutagenicity for Compound 4808 was 0.000 (discriminant score = −29.801). The computed rat oral lethal dose 50 (LD 50 ) value for Compound 4808 obtained with the Rat Oral LD 50 heteroaromatic model was 111.8 mg/kg. The lower and upper 95% confidence limits were 13.5 mg/kg and 925.8 mg/kg, respectively, and the computed rat oral LD 50 log value (1/Moles) was 3.630. These parameters were checked to be within their standard ranges and compound 4804 is suitable for further development as a lead compound.

Conclusions
Seasonal influenza is considered to be a serious public health problem. Drug resistance to existing influenza drugs requires researchers to develop new anti-influenza mechanisms. PB1-PB2 interactions are very important for the formation of RdRp, which is a key structural component of virus replication, so it is an attractive target for drug development. To date, no specific inhibitor of PB1-PB2 interactions has been reported. In this study, our main goal was to obtain a lead compound with drug-forming potential through a virtual screen. Firstly, according to the analysis of the Ramachandran diagram, protein 2ZTT is chosen as the docked protein structure. A pocket of PB1 protein, which currently has no drug development in this pocket, was selected as a target for drug design. In order to ensure a low false positive of the screened drug, we chose a compound library targeting PPI inhibitors and a cross-validation scheme using three software tools for docking. Finally, three compounds were screened out through molecular docking. In order to further ensure the effectiveness of screening, molecular dynamics simulation was performed for 10 ns. The simulations results show that these compounds could be stably bound in PB1 pocket, and the binding effect is more than 100 kJ/mol. We chose compound 4808 as the screening result, and analyzed the free energy decomposition of compound 4808. Based on the decomposition results, it is possible to improve drug activity by reducing the bad contact with protein residues 30, 41 and 45, and further enhancing the binding effect of residues 15, 59, 62 and 66. However, before entering clinical practice, further work still needs to be done on PB1-PB2 inhibitors. For example, there are few co-crystallization studies on PB1-PB2 inhibitors. The activity of them only reaches the micromolar level, and further structural optimization is required. This study provides a lead compound for developing new targets of antiviral drugs, so as to solve the drug resistance problem of current anti-influenza drugs, and puts forward suggestions for improving the compound. However, the compound was only computationally compared with existing active compounds, and further studies are needed to obtain in vivo activity test data.

Funding:
The authors gratefully acknowledge financial support from the Graduate Innovative Fund of the Graduate Innovative Fund of Wuhan Institute of Technology (No. CX2020022). The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.