Hydrogen bonding penalty used for virtual screening to discover potent inhibitors for Papain-Like cysteine proteases of SARS-CoV-2.

The Papain-Like proteases (PLpro) of SARS-CoV-2 play a crucial role in viral replication and the formation of nonstructural proteins. To find available inhibitors, the 3D structure of PLpro of SARS2 was obtained by homologous modelling, and we used this structure as a target to search for inhibitors through molecular docking and MM/GBSA binding free energy rescoring. A novel hydrogen bonding penalty was applied to the screening process, which meanwhile took desolvation into account. Finally, 61 compounds were acquired and 4 of them with IC50 at micromolar level tested in vitro enzyme activity assay, which includes clinical drugs tegaserod. Considering the importance of crystal water molecules, the 4 compounds were re-docked and considered bound waters in the active site as a part of PLpro. The binding modes of these 4 compounds were further explored with metadynamics simulations. The hits will provide a starting point for future key interactions identified and lead optimization targetting PLpro.

INTRODUCTION

There is a new public health crisis threatening the world with the emergence and spread of SARS‐CoV‐2 (SARS2) (Shin et al., 2020; Singhal, 2020). An early and essential process of SARS2 is the cleavage of polyprotein into 16 mature components termed nonstructural proteins (NSPs). Papain‐Like protease (PLpro) catalyzes itself and other NSPs release from the polyprotein, thereby initiating virus‐mediated RNA replication (Klemm et al., 2020; Ratia et al., 2008). The same as SARS‐CoV (SARS) PLpro, SARS2 PLpro has deubiquitinating and deISGylating enzyme activities, and evades the innate immune response of host cells by deubiquitinating and deISGylating structural and nonstructural proteins of SARS2 (Klemm et al., 2020; Mielech, Chen, et al., 2014; Mielech, Kilianski, et al., 2014). PLpro uses the thiol group of cysteine as a nucleophile to attack the carbonyl group of the scissile peptide bond (Chou et al., 2014). In the cytosol, the membrane‐associated PLpro domain found in the boundaries of nsp1/2, nsp2/3, and nsp3/4, recognizes the P4–P1 consensus cleavage sequence LXGG, and proteolytic of the peptide bond occurred after the glycine at the position, which is an essential process for viral replication (Gao et al., 2021; Han et al., 2005; Rut et al., 2020). The numerous functions and requisite roles of PLpro in viral replication and pathogenesis suggest that PLpro is essential for the viral life cycle, thereby it is a vital target for drug discovery studies against the recent epidemics caused by SARS2 (Klemm et al., 2020).

GRL0617 is a potent, selective, and competitive noncovalent inhibitor of SARS PLpro/deubiquitinase, with an IC50 of 0.6 μM and K i of 0.49 μM (Ratia et al., 2008). Aiming at the GRL0617 binding site, structure‐based inhibitors designs have been reported (Freitas et al., 2020; Fu et al., 2021; Rut et al., 2020). The main body of this site is populated with polar amino acids, with only one end consisting of two prolines forming a hydrophobic region (Freitas et al., 2020; Fu et al., 2021; Gao et al., 2021; Rut et al., 2020). Water‐mediated hydrogen bonds have been proposed to play important roles in PLpro–ligand complexes. In crystal structure 7JIT (C111S mutant), a water molecule mediated hydrogen bond between Lys157 and carbamylurea of Y95, a similar phenomenon is also seen in crystal structure 7JIW (wild‐type). Not only that, the sidechain of Glu167 showed minor migration due to water molecule extrusion, which resulted in Glu167 forming a new hydrogen bond with ligand (Shen et al., 2021). Water molecules in the binding pocket can not only occupy a certain space but also stabilize the residues around them by complexation (Osipiuk et al., 2021; Trujillo et al., 2012). Improper treatment of water molecules leads to a dramatic decrease in inhibitor activity, in contrast, the reasonable design of molecules to replace the water molecules in the binding pocket can improve the activity (Fell et al., 2020; Trujillo et al., 2012). The process which drugs exert their effects involves competition between ligands and water molecules in the binding pocket of the target protein. In other words, if a ligand wants to drive away water molecules and occupy the active pocket, it must break hydrogen bonds between water molecules and amino acid residues and form new hydrogen bonds (Abel et al., 2008). Hence, it is very important to correctly handle the hydrogen bond interaction between protein and ligand. However, in virtual screening methods, the scoring function stipulates the solvent as a continuous medium model that is poorly suited for bound water and therefore tends to give poor results in predicting binding affinities (Feig et al., 2004; Friesner et al., 2006; Schneider, 2010). To this end, there have also been attempts to consider bound water as part of the receptor for molecular docking, with statistical results showing improved precision (Thilagavathi & Mancera, 2010). Unfortunately, MMPB/GBSA does not deal well with water molecules in the receptor structure (Hayes et al., 2011), usually the water molecules in the crystal structure need to be removed for subsequent calculation, which reduced the accuracy of the scoring function.

To address these concerns, we applied the hydrogen bonding penalty (HBP) to the calculation of binding free energy (Zhao & Huang, 2011). Hydrogen bond penalty is based on a theory that hydrogen bonds formed between receptor and water molecules are considered optimal status, and after a water molecule is replaced by a ligand, energy loss will occur if the new hydrogen bonds energies cannot well compensate for the original hydrogen bonds formed between water molecular and receptor. Due to the large degree of freedom of water molecules, the entropy gain caused by the substitution of ligands is beneficial to the enhancement of the binding affinity between receptors and ligands. In contrast, usually, the enthalpy gain cannot be fully complemented because the newly formed hydrogen bonds do not conform to reasonable geometric constraints or cannot form new hydrogen bonds after desolvation, regardless of ligands or receptors (Mahmoud et al., 2020; Young et al., 2007). Therefore, we use HBP to characterize the enthalpic penalties. Lewater, a program for hydrogen bond penalty, is designed to predict the position of crystal water and polar interactions at the protein–ligand interface. The program calculates the enthalpy by detecting both the solvation and hydrogen bonding state at the protein–ligand interface. In recent years, hydrogen bond energies have been introduced into modified binding free energy calculations. Hydrogen bond energy usually consists of the distance function and angle function between hydrogen bond donors and hydrogen bond acceptors (Bao et al., 2020; Zhao & Huang, 2011). Solvent exposed moieties of ligand may also form hydrogen bonds with surrounding residues, but these interaction sites are competed by the solvent, consequently, estimating whether these moieties contribute to protein–ligand binding or given a weight value is necessary, which is commonly characterized by the ratio of solvent accessible surface area (SASA) (Mahmoud et al., 2020). Indeed, SASA is a very important parameter term when calculating hydrogen bond energy and hydrogen bond penalty because it directly affects the angle of hydrogen bonds and the binding “gesture” of ligands. Lewater takes into account both the polar and hydrophobic character of the binding site residues and ligands to predict the possible water molecules' binding points in the structure and makes a reasonable treatment for the sites where hydrogen bonds formed. Additionally, if the hydrogen bonds between receptors and donors are formed within solvent inaccessible regions, nonhydrogen bond penalty is created (Bao et al., 2020; Wang et al., 2002; Zhao & Huang, 2011).

METHODS

Sequence alignment and homologous modelling

Before the crystal structure of SARS2 PLpro was solved, we build the structural model using the homology modelling method. The sequence information between SARS2 and SARS, including secondary structure distribution, disulfide bond position, highly conserved residues, and partially conserved residues alignment were analysed by ESPript (Robert & Gouet, 2014). The 3D model was obtained by homologous modelling with SWISS‐MODEL (Waterhouse et al., 2018) using the crystal structure of PLpro of SARS (PDBID: 3E9S) as a template.

Screening process

The receptor model (model_01.pdb) and ligands were optimized by correcting the bond order, adding hydrogen atoms, distributing charges, and predicting the protonation states (pH 7.0). The ligand library is taken from “TargetMol, USA” (ChemDiv and “Bioactive Compound Library”). From 9175 compounds, 656 compounds with ligand efficiency (LE) <−0.3 and score <−8.0 kcal/mol were obtained by molecule docking using the professional version ledock (ledock_go) (Wang et al., 2016), then we separate the first pose of docking results of each ligand and apply prime‐mmgbsa (Jacobson et al., 2004) to each complex to calculate binding free energy (ΔG, kcal/mol). After integrating HBP into binding free energy calculation, a value of 1 in enthalpy calculated by Lewater corresponds to about 1.7 kcal/mol in binding free energy (Zhao & Huang, 2011) which is consistent with the experimental value (Fersht et al., 1985), 61 compounds were obtained as active molecules, the screening process can be seen in Figure 1.

FIGURE 1

Virtual screening workflow

Protein expression and purification and inhibitor activity detection

The catalytic domain of SARS2 PLpro used for in vitro compound adduct formation study was sub‐cloned into pET‐28a vector. This bacterial‐expression construct encoded a fusion protein with N‐terminal 6 × His, protein expression was performed in Escherichia coli BL21 (DE3) by induction with 0.1 mM IPTG at 18°C overnight, and purification was achieved using Ni‐NTA chromatography. Cells were harvested after 16 h of incubation at 20°C and lysed by sonication in 20 mM phosphate buffer containing 500 mM NaCl, 200 μg/ml of lysozyme, and 0.1 mM PMSF, pH 7.4. Cell lysates were clarified by centrifugation at 13,201 g, 4°C, for 15 min and applied to a Ni‐NTA resin. Recombinant protein PLpro was eluted with binding buffer containing 5 mM imidazole and extraction buffer containing 250 mM imidazole, respectively. PLpro was collected in 4 fractions, 5 ml each. The OD280 of the 4 fractions were detected by NanoDrop Spectrophotometer and the fractions with concentration below 0.1 mg/ml were discharged. To further concentrate the PLpro, the protein elution was transferred to a Millipore Amicon Ultra‐15 centrifugal filter tube immediately to centrifuge at 4500 g at 4°C until there was 5 ml protein left. The protein was diluted with assay buffer containing 200 mM NaH2PO4, 200 mM Na2HPO4, and pH 6.8.

The purified protease ran at around 36.5 kDa on sodium dodecyl sulfate‐polyacrylamide gel electrophoresis, with over 90% purity (data not shown). To find an effective inhibitor, 200 nM PLpro and 20 μM fluorogenic substrate (MCA‐Dnp‐Lys) in buffer (20 mM Na2HPO4–NaH2PO4 at pH 6.8) were mixed, and different concentrations of the compounds were added at the incubation temperature of 25°C, shock incubation at room temperature for 15 min, then quickly add 20 μmol/L fluorogenic substrate1 (MCA‐Dnp‐Lys) which sequence was derived from the amino acid sequence of the viral structural fragment NSP2/3 identified by PLpro. The enhanced fluorescence emission upon substrate cleavage was monitored at the excitation and emission wavelengths of 320 and 425 nm, using SpectraMax iD5, respectively, in the presence or absence of the compounds assayed. Because the active site of PLpro is similar to that of ubiquitin, Cbz‐RLRGG‐AMC which sequence is derived from five c‐terminal residues of ubiquitin was used as the fluorogenic substrate2, 40 nM PLpro and 2.4 μM fluorogenic substrate in buffer (20 mM Na2HPO4–NaH2PO4 at pH 6.8) were mixed, different concentrations of the compounds were added at the incubation temperature of 25°C, shock incubation at room temperature for 15 min, then quickly add 20 μmol/L fluorogenic substrate2, the enhanced fluorescence emission upon substrate cleavage was monitored at the excitation and emission wavelengths of 360 and 460 nm, using SpectraMax iD5, respectively, in the presence or absence of the compounds assayed. The logarithm of compound concentration was taken as the abscissa, and the corresponding inhibitory rate was taken as the ordinate. The IC50 value of the compound was calculated by the four‐parameter method.

Docking with explicit water molecules

To evaluate the role of ordered water molecules, the crystal water molecules in the binding pocket were included in the docking. We modelled the ordered water molecules based on the template X‐ray structure (PDB ID: 3E9S). Five crystal water molecules were included and hydrogen atoms were generated with Chimera (Pettersen et al., 2004). All other parameters were consistent with the docking method in the screening process. The top 3 docked poses for each compound were obtained for further analysis.

Conformational verification

To hunt for the correct binding poses, we used metadynamics (Cutrona et al., 2020) to evaluate the docking results. The receptor structure containing 5 crystal water molecules was prepared using Protein Preparation Wizard panel (Sastry et al., 2013), including adding hydrogen atoms, removing heteroatoms, and protonation state prediction (pH 7.0). The docked poses were evaluated with the Binding Pose Metadynamics. The stability of the binding pose was evaluated with the RMSD over the course of the simulation and the persistence of key interactions between the ligand and the receptor. The docked poses from ledock_go for each ligand were evaluated with 10 trials each. To facilitate differentiation, we created an input file for each ligand.

Molecular dynamics simulation

We employed the Desmond (Maestro‐Desmond Interoperability Tools, Schrödinger, New York, NY, 2020‐2) program to run molecular dynamics simulation based on the OPLS3e force filed (Harder et al., 2016). The protonation states of residues in PLpro and ligands were assigned at pH 7.0. The complexes were solvated in a cubic box of TIP3P water models and NaCl was added to neutralize the system. The system was subjected to 100 ps minimization (10,000 steps with Brownian motion simulation) to adequate equilibrated complexes and solvent molecules. The harmonic position restraints were applied on the backbone of PLpro with a force constant of 5 kcal.mol−1.Å−2. A 5 ns constant volume and constant temperature (NVT ensemble) simulation at 300 K was carried out to allow adequate equilibration of the water molecules with harmonic position restraints were applied on the backbone (force constant 10 kcal.mol−1.Å−2) and side chains (force constant 5 kcal.mol−1.Å−2). The 20 ns constant pressure and constant temperature (NPT ensemble) simulations at 300 K and 1.01325 bar were carried out with 5 kcal.mol−1.Å−2 force constant restraint for backbone, and 100 ns equilibrium MD simulation was performed with NPT ensemble at 300 K with a Nose–Hoover chain thermostat and 1 atm with an isotropic Martyna–Tobias–Klein barostat. The time step was set to 2.0 fs. Short‐range electrostatics was calculated with a cutoff of 9.0 Å and long‐range electrostatics was calculated with QuadS. The trajectories were saved every 20 ps for further analyses.

RESULTS AND DISCUSSION

Structural comparison

PLpro of SARS2 and SARS are highly homologous (their identity is 82.9% and similarity is 94.9%) Their active sites are very similar. The details can be seen in Figure 2. The X‐ray structures of SARS and SARS2, and the SARS2‐PLpro homology model superimpose very well (with binding site atom RMSD 0.72 Å). The key residues in the active pocket adapt highly similar structure except for Glu167 and Leu162 (Figure 3). It is known that the side chains of Glu167 and Leu162 are highly flexible, were there to fit is typically affected by ligand (Figure S1 and Table S1) (Gavory et al., 2018; Leger et al., 2020; Wang et al., 2020).

FIGURE 2

PLpro's sequence alignment results of SARS‐CoV‐2 (SARS2) and SARS‐CoV (SARS). α, β, and η are symbols of the secondary structure of SARS2, the arrows represent β‐strands, the α represents α‐helices, the η represents 310 helices, and “TT” represents strict beta turns. White characters with a red background are the completely conserved amino acids. The red character with white background means the consensus is >70, the black character with white background means the consensus <70. Block uppercase is identity, “!” is anyone of IV, “$” is anyone of LM, “%” is anyone of FY, and “#” is anyone of NDQEBZ. Blue dots represent residues in the active pocket of SARS2 and SARS

FIGURE 3

Homologous modelling result. Model‐01 is the structure model, 3E9S is the crystal structure of SARS and 7CJM is the crystal structure of SARS2. N is the N‐terminal and C is the C‐terminal. The red box represents the binding site region, and the snapshot on the bottom right shows the distribution of residues in the active pocket

Virtual screening and enzyme activity experiments

The potential binding sites of water molecules in the active pocket were predicted with Lewater (Zhao & Huang, 2011) (for compound 56 shown in Figure 4). Obviously, as the ligands occupy positions of four crystal water, water molecules must be discharged into the solvent, resulting in energy loss. For polar residues around the ligand, all the potential points form a surface, which indicates bound waters distributed at these areas. These regions are exposed to the aqueous solution, and there would be an energy penalty in hydrogen bonds with the ligands if the hydrogen bonds were to exit. Therefore, after binding free energy to rescore by MMGBSA, HBP was used to correct the score values. From 61 hit compounds, 4 of them shown inhibitory activity. The results can be seen in Table 1 and Figures S2 and S3.

FIGURE 4

The area of bound water is predicted by Lewater. The red surface represents the possible location of water molecules. The carbon atoms of the ligand and protein are green and grey, the oxygen atoms are red, and the nitrogen atoms are blue, crystal water molecules were yellow spheres. The ligand is compound 56

TABLE 1

Chemical structures and inhibitory activity of the SARS‐CoV‐2 PLpro inhibitors. MCA‐Dnp‐Lys and Cbz‐RLRGG‐AMC are different substrates for SARS‐CoV‐2 PLpro

Compounds	Structure	Score	LE	MMGBSA	HBP	IC50 (μM)	IC50 (μM)
25		−9.0	−0.30	−65.6	2.7	45.90	60.28
27		−9.5	−0.31	−59.9	2.4	31.03	59.09
38		−8.1	−0.32	−59.4	2.0	19.67	>100
56		−8.7	−0.39	−62.1	3.1	9.25	1.42

Note: Score was calculated by ledock_go, the unit of Score, MMGBSA, and HBP is kcal mol−1. LE is the abbreviation of Ligand efficiency. Column 7 is the IC50 measured with MCA DNP Lys as the substrate and Column 8 is the IC50 measured with Cbz‐RLRGG‐AMC as the substrate.

Poses evaluation

We evaluated the effects of explicit water molecules in docking. With the crystal water molecules included, the docking poses were closer to that of the GRL0617 in the crystal structures 3E9S and 7CJM (Figure 6a–d). Superimposing the poses of 4 compounds with GRL0617 (Figure 6e), a similar trend can be found—they form a “V‐shaped” structure, and the turning point occurs at “knee point nitrogen atoms” (Figure 6e, kermesinus oval frame) in the middle of the straight‐chain (compound 25 locate at piperazine ring). Interestingly, the coordinates of this knee point are very close, which indicates that the nitrogen atom here is very important for the inhibitory activity of the inhibitor (Fu et al., 2021; Gao et al., 2021).

FIGURE 6

(a) Compound 25. (b) Compound 27. (c) Compound 38. (d) Compound 56. Superimposing poses of regarded crystal water molecules as a part of the protein or not with GRL0617 in 3E9S and 7CJM. (e) Superimposing re‐docking poses of 4 compounds with GRL0617 in 3E9S and 7CJM

Molecular dynamics simulation analysis

The binding model indicates that the cavity of PLpro is filled with polar residues (Figures 7 and 8). The compound 25 binding with PLpro by π‐cation with residue Tyr264. Compound 27 forming salt bridge (with residue Asp164), water bridge (with residue Glu161), π‐cation (with residues Tyr264 and Lys157), salt bridge (with residue Asp164), and hydrogen bonds (with residues Asp164 and Gln269) with PLpro. There are hydrophobic interactions (with residues Tyr264 and Pro248, not displayed), hydrogen bonds (with residues Gln269), and Pi–Pi stacking (with residue Tyr268) between compound 38 and PLpro. As for compound 56, salt bridge (with residue Asp164), π‐cation (with residue Tyr264), water bridge (with residue Ala246 and Glu167), and hydrogen bonds (with residues Asp164, Gln269, and Tyr273) in the active pocket of PLpro. From all residues marked in Figures 7 and 8, Tyr264 had the highest frequency between 4 compounds and PLpro, followed by Gln269 (compound 27, compound 38, and compound 56) and Asp164 (compound 27 and compound 56), which suggests that they have a relatively important role in maintaining affinity.

FIGURE 7

Interactions between 4 compounds and PLpro analysed by MD (more than 30% rate is displayed). The percentage represents the probability of total MD process. (a) compound 25, (b) compound 27, (c) compound 38, and (d) compound 56

FIGURE 8

3D modes between 4 compounds and PLpro. (a) compound 25, (b) compound 27, (c) compound 38, and (d) compound 56

In addition to protein‐ligand interactions, the internal energy of the ligand is an important factor in assessing binding affinity (Kuhn et al., 2016; Pennington & Moustakas, 2017; Senger et al., 2007). Torsion energy analysis shows that these compounds have 7, 9, 6, and 9 rotational bonds, respectively (Figure S4). The diphenyl ether of compound 25 (Figure S4A) in an extremely flexible state, with the dihedral angle (red and green) sampling extensively between −180° and 180°, the chlorobenzene (deep purple), diethylamine (cyan and pink), and aliphatic chain (wheat, purple, and lime green) of compound 27 (Figure S4B) and aliphatic chain (pink, cyan, and lime green) and methoxy (deep purple) of compound 56 (Figure S4D) sampled much low energy poses, which result in higher RMSD (Figure S5) than compound 38 (Figure S4C). Rotatable bond, linking the “knee point nitrogen atoms” with others, centrally distributed at high energy dihedral angle states for compound 25 (lime green, about 3.68 kcal/mol) and compound 56 (purple, about 4.40 kcal/mol).

The rotatable bonds of compound 25 revealed high entropy (red and green, Figure S4A) and torsional potential energy (lime green, Figure S4A) after entering the binding site. Trajectory analysis shows that compound 38 and compound 56 have lower RMSD (Figure S5) and RMSF (Figure S6). During the simulations, compound 25 caused large fluctuations in the protein backbone (Figures S7 and S8), compounds 27 and 56 had comparable effects, and compound 38 caused the smallest fluctuation.

In the statistics of metadynamics interactions, compound 25 had many interactions with PLpro (Table S2), however, an unstable binding pose leads to a low proportion of many interactions in a long time MD process. For example, in the first 10 ns (0–500 frames) of the MD process, the hydrogen bond formed with Gln269 appeared 365 times, accounting for 73.00%, but 1334 times in 100 ns (5000 frames), accounting for 26.68% (data not shown), which further explains the reason why the low activity of compound 25. As for compound 27, π‐cation interactions with Lys157 appeared 223 times in the first 500 frames, and 4198 times in the whole MD process (5000 frames). Similarly, Tyr264 appeared 369 times and 4183 times in the first 10 ns and all simulation time (100 ns), respectively (data not shown), which indicates that the contribution of these interactions was relatively weak in the early MD but their importance were reflected along with the simulation system gradually stabilized. The same reason can be invoked to explain the case of hydrogen bonds between compound 56 and Asp164 and Gln269 (data not shown), as the results of metadynamics simulations (Table S2), the average value of these interactions in the 10 evaluation results was very low. Although compound 27 and compound 56 were comparable in terms of interaction, the RMSD (Figure S5) and RMSF (Figure S6) of compound 27 were higher than compound 56, and metadynamics analysis (Figure 5 and Table 2) also confirmed, in binding poses, compound 56 (56‐3) being more stable than compound 27 (27‐2). On the contrary, compound 38 showed the lowest RMSD and RMSF, however, a paucity of interactions resulted in its affinity being much lower than that of compound 56. It is worth mentioning that the torsional potential energy of ligands is a potential factor to improve binding affinity (Engelhardt et al., 2019; Muegge et al., 2017; Sellers et al., 2019) and the potence of compound 56 will be better if the dihedral angle of rotatable bond (purple, Figure S4D) been optimized.

FIGURE 5

Comparison of metadynamics results among 4 compounds. The abscissa is the simulation time of molecular dynamics, and the ordinate is the mean value of RMSD in 10 trials of the compound

TABLE 2

Metadynamics analysis of docking poses for 4 compounds. The persistence of Pi–Pi stacking was zero (not shown)

Number	PoseScore	Persistence	Persistence HBond	CompScore
25‐1	2.81	0.36	0.36	1.01
25‐2	2.99	0.57	0.57	0.13
25‐3	3.21	0.53	0.53	0.55
27‐1	2.91	0.36	0.36	1.12
27‐2	2.38	0.36	0.36	0.56
27‐3	4.41	0	0	4.41
38‐1	6.55	0.05	0.05	6.30
38‐2	1.61	0.52	0.52	−0.99
56‐1	2.79	0.08	0.08	2.38
56‐2	2.65	0.24	0.24	1.47
56‐3	2.28	0.63	0.63	−0.86

Note: PoseScore was the expectation of the RMSD of the pose over the course of the metadynamics. Persistence was the average persistence of contacts over the course of the metadynamics trajectories. Persistence HBond was the average persistence of HBonds. CompScore was linearly combining the PoseScore and Persistence scores.

Tegaserod, a 5‐hydroxytryptamine (serotonin) type 4 receptor agonist, is used to treat lower gastrointestinal hypomotility associated with irritable bowel syndrome (IBS‐C), and exhibits favourable efficacy for the treatments of colitis (Madia et al., 2020). From Figure 7d, the precursor of pentylaminoguanidine (PAG) undertakes the main interactions between tegaserod and PLpro. In previous research, PAG, a product of acid hydrolysis of tegaserod in the stomach, is the crucial consistency that exerts pharmacodynamic effects for IBS‐C. SARS2 can cause a severe inflammatory response, which is characterized by the continuous increase of inflammatory cytokines, such as IL‐6 and TNF‐α (Del Valle et al., 2020). IL‐6 can activate JAKs (Janus kinase), which will phosphorylate signal transducer and activator of transcription 3 (STAT3), a crucial transcriptional factor that is overactivated in multiple cancer types and many other diseases. However, tegaserod can down‐regulate STAT3 activity by specifically inhibiting the phosphorylation of STAT3 in cancer cells, and this inhibition is independent of the 5‐HT4R (Zhang et al., 2020). Thus, tegaserod has potential activity against SARS2 but also has anti‐tumor activity, which is good news for cancer patients with SARS2 infection. What's more, inhibition of PLpro activity also contributes to the innate immune function of the host, which is dual effects on the treatment of viral infection (Klemm et al., 2020; Mielech, Chen, et al., 2014; Mielech, Kilianski, et al., 2014). In summary, tegaserod showed a prominent inhibitory effect on PLpro. As a clinical drug, tegaserod may be a good choice for anti‐COVID‐19.

CONCLUSION

In this work, virtual screening and molecular dynamics simulations were used to discover novel inhibitors for PLpro. Among 61 compounds tested, 4 of them were found to have micromolar inhibitory activity against PLpro. Molecular docking results showed that the existence of bound water molecules had crucial effects on the binding conformation of ligands, and the “knee point nitrogen atoms” is the distinct features of active ligands. The structure–activity relationship suggests that Tyr264 and Gln269 are key residues for PLpro.

CONFLICT OF INTEREST

There is no conflict of interest.

Figures S1‐S8

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.