2.1 Virtual Screening
Protein Activity Site Analysis
In this study, we obtained the Mpro structure from the PDB protein database (https://www.rcsb.org, ID: 6LU7), and the native ligand for Mpro was the potential molecular inhibitor N3 [18]. N3 is a broad-spectrum inhibitor targeting the main protein of coronavirus designed by Haitao Yang et al. in 2005 [25]. By binding to the active site of the SARS-CoV-2 main protease and interacting with the Cys145 residue through Michael addition reaction, N3 is capable of forming a stable covalent bond that exerts an inhibitory effect on the activity of the protein. The three-dimensional spatial structure and secondary structure of the main protein can be observed in 6LU7 (Figure 2A, B) [26-29], which consists of two domains, including 10 α-helices and 13 β-sheets. The ligand binds to the active site near β3, β4, and β16 with Cys145 and His41 as key residues.
SP Docking
We conducted standard precision molecular docking (SP docking) on all 320 molecules, including the positive compound Nilotinib, to assess their affinity to Mpro (Table S1) within active site. The molecular structures of the compounds were drawn in 2D Sketcher of Maestro (version v128117, Schrödinger 2018-1, New York, USA) and subsequently prepared for docking through the use of LigPrep (LigPrep, Schrodinger, LLC, New York, NY, USA). Prior to docking, the protein underwent preparation and energy optimization with the aid of Wizard in Glid. The receptor grid generation tool was employed to generate a grid box centered on native ligand N3, with adjustments made to the grid size to optimize calculations. The minimized molecules were then docked to the grid box, utilizing standard precision docking. The results indicated that most drugs had binding energy between -5.5 to -7.5 kcal/mol (Figure 2), while 63 compounds exhibited higher docking binding energy than Nilotinib, with ΔG > -7.026 kcal/mol. These values corresponded to their docking binding affinity (Kd) > 142303.6 M-1 (Table 1 and S1).
<Figure 2>
XP Docking and MM-GBSA Evaluation
As extra precision docking (XP docking) mode is more advanced, it helps filter out false positives and provides a better association between docking score and good poses [30]. Therefore, XP docking was employed to re-dock the top 64 molecules, including nilotinib, with Mpro and rank compounds by their XP binding energy (Table 1). Based on the protein and re-docked ligands, Molecular Mechanics/Generalized Born Surface Area (MMGBSA) analysis was conducted to evaluate the free ligand binding energy of 48 compounds with higher XP binding energy than Nilotinib (ΔG > 5.476 kcal/mol). The results of analysis identified 32 molecules with higher free binding energies compared to Nilotinib (dG > -48.96 kcal/mol), which have been selected for further screening studies.
<Table 1>
2.2 Further Screening through Infuced Fit Docking
The induced-fit theory was introduced by Koshland in 1995, which posits that enzymes do not have a complementary shape to the substrate but form a complementary shape only after induction [31]. When it comes to ligand-receptor interaction, conformation of the receptor, especially around the binding site, was also induced to be altered, better matching the shape and binding pattern of ligand molecule. In this study, to obtain more realistic models of molecule-protein interactions, Induced Fit Docking (IFD) is applied for further virtual screening.
The 32 compounds were chosen to re-docked into Mpro by applying Induced Fit Docking (IFD), with Nilotinib as control, and calculated their Induced Fit Docking binding energy. Multiple docking results are produced during IFD docking due to the various poses of the compounds, and we recorded the results with the highest docking energy of each compounds in Table 1. The interaction of protein-ligand complexes were then displayed and analyzed in Maestro and potential interaction of all 32 compounds were counted (Figure 3 – 8 and Table S2). In comparison to SP or XP docking, IFD docking facilitated stronger binding and more protein interactions with filtered molecules, further accentuating differences in affinity between these compounds and Mpro. As shown in Table 2 (or Table S2), Nilotinib interact with Mpro mainly through hydrogen bonds, consistent with that reported in literature. For the screened compounds, they are able to engage more amino acids in protein-ligand interactions and display new interactions such as electrostatic interaction, pi-pi stacking hydrophobic interaction, compared with Nilotinib.
Here, considering the diversity of structural features and binding interaction models, 5 molecules are selected as representatives to illustrate the molecule-protein interactions. Nilotinib was employed as a control to further clarify the mechanism behind the compounds' capacity to generate stronger interaction forces.
<Table 2>
Nilotinib and Mpro Interaction
Nilotinib can bind to Mpro well in its active pocket (Figure 3). The complex consists of seven hydrogen bonds (Glu142, Cys145, His164, Glu166, Gln 189) including 1) two N atoms on the pyrimidine with Cys145, Asn142 and Gln189, respectively; 2) NH linked to the pyrimidine with Asn142; 3) CO in amide with Asn142 and Gln189; and 4) NH in amide with Glu166. The IFD docking binding energy of stabilization is ΔG = -9.179 kcal/mol.
<Figure 3>
Compound V253 and Mpro Interaction
The compound V253 has a pyridine biphenyl amide structure with serine as Linker. It could form ten hydrogen bonds (Leu141, Gly143, Ser144, Cys145, His164, Glu166, Gln189 and Gln192), two pi-pi stacking (His41) and a salt bridge (Glu166), which enhanced the stability of the complex. In addition to the similar hydrogen bonds as Nilotinib, there are also other hydrogen bonds, including: 1) the NH of amide in serine backbone with His164 and Gln189; 2) the OH in serine side chain with Leu141; 3) the N atom on pyridine with Gln192; and 4) the NH of amide linked to the pyridine with Glu166. Notably, V253 is capable of forming two pi-pi stacking interactions with the imidazole ring of His41. Furthermore, the acyl ethylenediamine connected to the pyridine contains a primary amine which have been protonated, allowing it to form salt-bridge interactions with Glu166. Overall, V253 could bind to Mpro through hydrogen bond and pi-pi bond with His41 and Cys145, thereby impairing the activity of protease.
<Figure 4>
Compound V247 and Mpro Interaction
The structure of Compound V247 closely resembles that of V253, with the only difference being the presence of a tBu protecting group on the serine sidechain of V247, thus making the docking pose and ligand-protein interactions similar to that of V253. The docking result of compound V247 with Mpro showed six hydrogen bonds (Asn142, His164, Glu166, Gln189 and Gln192) and a salt bridge with Glu166 due to protonated ammonia in ligand, thus forming a stable complex. Exclusively, the CO of amide in V247 generate one hydrogen bond with Asn142, which was different from that of V253. Interestingly, despite having fewer interactions, V247 demonstrates similar binding energy compared to V253.
<Figure 5>
Compound V133 and Mpro Interaction
The V133 shares a comparable structure with compound V253, except for the linker type of amino acid present in its molecule, which is alanine instead of serine. Thus, the interactions type between V133 and Mpro was quite similar to that of V253. As shown in the figure 6, compound V133 binds to the protein through six hydrogen bonds (Gly143, Cys145, Glu166, Arg188, and Gln189), and further promotes the stability of the complex through the salt bridge formed between protonated ammonia (primary nitrogen in ethylenediamine) in ligand and Glu166. Uniquely, the CO of amide attached to the pyridine demonstrates two hydrogen bonds with Gly143 and Cys145, respectively.
<Figure 6>
Compound V228 and Mpro Interaction
Compound V228 has backbone structure with tertiary leucine as Linker. The docking result shows that compound V228 can bind to the active site of Mpro, and their interaction includes seven hydrogen bonds (His41, Asn142, His164, Glu166, His172 and Gln189). Compare to V253, incorporation of tert-butyl into compound V228 confers significant steric hindrance, leading to a distinct conformation in the docking results. The amide structure attached to pyridine exhibits hydrogen bonding interactions with His41, His164, and Gln189, while another amide structure linked to the halogenated benzene engages in hydrogen bonding with Glu166 and Asn142. Unlike the above compounds, V228 with a cyclopropanamide structure lacks primary amine so that it cannot form a distinctive salt-bridging interaction with Glu166.
<Figure 7>
Compound V291 and Mpro Interaction
Notably, the pyrrolidine of proline in V291 provides the basis for restricting the deformation of molecule during binding to the protein, making its stable insertion and binding to the protein active sit. The indazole structure in V291 exhibits stronger hydrophobic interaction compared to the pyridine-linked pyrimidine structure in nilotinib, which facilitates ligand-protein interactions. This compound binds to Mpro through nine hydrogen bonds (Thr26, His41, Asn142, Gly143, Cys145 and Glu166), three Pi-pi staking (His41 and His163) and one salt bridge (Glu166). Indazole has the unique ability to form two pi-pi stacking interactions with His41, which is not found in compounds with other structural types. Hydrogen bond between ligand and Cys145 or His41 could prevent key residues cysteine and histidine from participating in substrate hydrolysis catalysis. Three pi-pi stacking improved the hydrophobic interaction to enhance the binding of complex.
<Figure 8>
Based on above results, it was found that these screened compounds interact with Mpro protein mainly through hydrogen bonding, pi-pi stacking, and salt bridges. Hydrogen bonding is the most common interaction in these molecule-protein complexes, and Gly143, Cys145, Glu166, and Gln189 are the key residues involved in hydrogen bonding interaction. In consistent with Nilotinib, some compounds, such as V133, V253 and V291, can directly interact with cysteine 145 through hydrogen bonding, leading to direct inhibition of the cysteine-mediated hydrolysis reaction. Besides, His41 is commonly involved in pi-pi stacking with aromatic structures in docking molecules. In particular, compound V291, with a larger indazole aromatic ring structure, exhibits stronger pi-pi stacking with His41 and hydrophobic effects, enabling His41 to participate more effectively in ligand-receptor interaction. Additionally, the secondary amine embedded in proline of V291 and the primary amines of other compounds facilitate the involvement of Glu166 in salt bridge generation between molecules and Mpro protein. Notably, for Cys145 and His41, two key amino acids involved in the enzymatic hydrolysis, Nilotinib only forms hydrogen bonds with cysteine, while the screened compounds can simultaneously interact with another key amino acid, His41, in pi-pi stacking, which might enhance enzymatic inhibition. Moreover, the reported compounds could generate salt bridge with Glu166, revealing a new bind mode with Mpro. These findings could provide novel perspectives for the development of Mpro inhibitors.
2.3 Pharmacophore Analysis
On the basis of molecular docking, we counted and analyzed the molecules which could bind with Mpro in the active site. Next, these molecules were selected for superimposition with applying their docking conformation (as shown in Figure 9A, B), and we further analyzed the pharmacophore of all molecules. As shown in Figure 9C, there are seven common components identified from superimposition: three Aromatic Rings (R10 dint R11, R12, orange), three Acceptors (A2, A3, A4, pink), and one Donor (D6, blue). Particularly, R10 and R11, A3, and D6 are deeply embedded in the pocket of the active site, thereby providing steric hindrance to substrate binding to the protein. Besides, A2 and A4 are located near the key cysteine amino acid to prevent the substrate from reacting with Cys145. Moreover, R12, which are frequently occupied by benzene rings, could exhibit hydrophobic interactions with the protein on the left side of the pocket. Meanwhile, R10 are predominantly captured by aromatic rings, allowing for effective pi-pi interactions with another crucial amino acid His41. In addition, the amide structures are mostly observed at A3 and D6, with electron-absorbing groups promoting the polarization of NH bond to be effectively acted as hydrogen bond donor, while the carbonyl oxygen acts as hydrogen bond acceptor, resulting in interactions with neighboring amino acids like Glu166. These findings are in agreement with the conclusions drawn from the ligand-receptor interactions demonstrated by the IFD docking results.
<Figure 9>
2.4 Molecular Dynamics Analysis
Virtual screening and molecular docking have shown the binding ability of different compounds to the main protease, while molecular dynamics simulations can provide insights into the dynamic stability of receptor-ligand complexes under physiological conditions. Thus, based on the binding affinity, ligand-receptor complex interaction and structural diversity of compounds (Figure S1), 20 compounds were selected to be further carried out molecular dynamics simulations to investigate their binding stability (Figure 10). Throughout the molecular dynamics simulation, one frame of trajectory was generated every 100 ps, resulting in 1000 frames after 100 ns of operation. Additionally, a thorough analysis of the resulting data was conducted, including Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), and Protein-Ligand Contacts analysis.
The conformation of the ligand-protein complex at 0 ns was used as reference for the calculation of RMSD values for all 1000 frames. During the simulation, 6 compounds (V131, V133, V172, V245, V247and V253) showed poor dynamic stabilities (Figure 10), while other compounds exhibited stable binding with Mpro protein in the active site. Fluctuations in RMSD values are usually associated with changes in ligand-protein interaction. For example, compound V247 displayed large variation of RMSD, which initially formed a strong bond with the Glu166 of Mpro, but weakened or even disappeared after 27ns (Figure S2). Figure S2 A shows the three-dimensional structure of the V247-Mpro complex at 1ns and 28ns, respectively. It is evident that the N-(2-(dimethylamino)ethyl)nicotinamide of compound V247 have detached from the protein, and the ethylenediamine portion has been entirely separated from the protein surface at 28ns, rendering it unable to interact with residues, thereby leading to strong fluctuations.
As for other 14 ligand-protein complexes with good dynamic stability, they showed stable RMSD values over long periods of time. To delve deeper into the dynamic changes in protein-ligand contact during the simulation, an analysis of Root Mean Square Fluctuation (RMSF) and Protein-Ligand Contacts was conducted (FigureS3 – S9). RMSF is useful for characterizing local changes along the protein chain, which helped identify the residues responsible for altering the fluctuations of the protein-ligand complex structure. For all complexes except the V222-Mpro complex and the V254-Mpro complex, the Protein RMSF values exhibited negligible variation, suggesting that the system was in a state of equilibrium throughout the simulation. The residues in contact with ligand were annotated in green, and the status of their contacts with the ligand was monitored every 0.2 ns throughout the 20 ns dynamic simulation. Most compounds maintain contact with a number of particular residues during the dynamic simulation, except for V75-Mpro complex, V97-Mpro complex and V282-Mpro complex.
<Figure 10>
2.5 ADME/T Prediction
ADME/T is a very important evaluation standard in contemporary drug design and drug screening. Any compound that exhibits drug-likeness must have moderate ADME/T properties. QIKPROP was used here to predict the ADME/T properties of the 20 compounds mentioned above (Table 3). The QIKPROP analysis includes the following standard limits: PSA (Van der Waals surface area of polar nitrogen and oxygen atoms and carbonyl carbon atoms), QPlogS (aqueous solubility, log S. S in mol dm-3 is the concentration of the solute in a saturated solution that is in equilibrium with the crystalline solid), QPlogPo/w (octanol/water partition coefficient), donorHB, accptHB , CNS (central nervous system activity), #metab (number of likely metabolic reactions), HumanOralAbsorption, QPlogBB (brain/blood partition coefficient), QPPMDCK (apparent MDCK cell permeability in nm/sec), QPPCaco (apparent Caco-2 cell permeability in nm/sec) and QPlogHERG (IC50 value for blockage of HEGR K+ channels). Based on the collective data, we infer that the ADME/T properties of ccompounds are within the prescribed limits for a potential candidate drug compound.
<Table 3>
2.6 Compound Enzymatic Activity Assay
Based on MD stimulations and ADME/T predicted results, 9 compounds are selected to assess the inhibitory potency and calculate IC50 values (Table 4). We used fluorescence resonance energy transfer (FRET)-based assay to measure these 9 compounds inhibitory activity against Mpro protein in virto. As shown in table 3, compound V291 exhibited a potent inhibitory activity with IC50 value of 2.77±1.56 μM, while other compounds did not exhibit significant inhibitory activity against Mpro within a concentration of 20 µM. This preliminary result provides biological evidence that compound V291, which has a similar structure to Nilotinib, could potentially exert an inhibitory effect on Mpro.
<Table 4>