Multi-conformation representation of Mpro identifies promising candidates for drug repurposing against COVID-19.

The COVID-19 main protease (Mpro), one of the conserved proteins of the novel coronavirus is crucial for its replication and so is a very lucrative drug target. Till now, there is no drug molecule that has been convincingly identified as the inhibitor of the function of this protein. The current pandemic situation demands a shortcut to quickly reach to a lead compound or a drug, which may not be the best but might serve as an interim solution at least. Following this notion, the present investigation uses virtual screening to find a molecule which is alraedy approved as a drug for some other disease but could be repurposed to inhibit Mpro. The potential of the present method of work to identify such a molecule, which otherwise would have been missed out, lies in the fact that instead of just using the crystallographically identified conformation of the receptor's ligand binding pocket, molecular dynamics generated ensemble of conformations has been used. It implicitly included the possibilities of "induced-fit" and/or "population shift" mechanisms of ligand fitting. As a result, the investigation has not only identified antiviral drugs like ribavirin, ritonavir, etc., but it has also captured a wide variety of drugs for various other diseases like amrubicin, cangrelor, desmopressin, diosmin, etc. as the potent possibilities. Some of these ligands are versatile to form stable interactions with various different conformations of the receptor and therefore have been statistically surfaced in the investigation. Overall the investigation offers a wide range of compounds for further testing to confirm their scopes of applications to combat the COVID-19 pandemic.

Introduction

COVID-19, the novel coronavirus (SARS-CoV-2), a member of the class β-coronavirus, is the causal agent of the viral pneumonia outbreak 2019–2020. The recent outbreak of this highly transmissible airborne disease, first officially recorded in Wuhan, China, and thereafter in many other countries all around the world, has been declared as pandemic by WHO in March 2020. As of January 13, 2021, more than 90 million people are reported to be affected worldwide [1], and one can suspect that the actual number of affected people is more than that as many asymptomatic, yet infected people are not yet tested. As the aforesaid virus is a novel member of the family Coronaviridae, until now, targeted therapies against the virus are scanty.

A number of proteins are encoded by the genome of this RNA virus including 3-chymotrypsin-like protease (3CLpro) or the main protease (Mpro), RNA-dependent RNA polymerase, helicase, etc. [2]. Among these proteins, Mpro becomes a very potential drug target because of a number of reasons [3, 4]. Replication and transcription of the virus are mediated by two overlapping polyproteins pp1ab. A highly complex proteolytic processing event of the polyproteins is catalyzed by two proteases. Among the two proteases, Mpro is the key enzyme which catalyzes the cleavage of the polyproteins at 11 sites conserved in SARS-CoV and SARS-CoV2. Therefore, Mpro is essential for viral replication and is a lucrative drug target [5]. Another important feature is that though the virus is very much prone to mutations, Mpro, being a conserved protein throughout the Coronavirinae subfamily is much less susceptible to frequent mutation. Therefore, designing inhibitors against Mpro would be worthwhile for the treatment of not only the SARS-CoV-2 but also for the entire subfamily. Also, cleavage specificity of no human protease is known to match with that of Mpro, thereby low or no inhibitor toxicity for human proteases can be expected. The Mpro consists of 3 distinct structural domains made up of 13 β-sheets and 11 α-helices (Fig. 1) [6]. Domain1 (residues 8–101 shown in orange) and domain 2 (residues 102–184 shown in green) comprised mostly of antiparallel β-barrel, with a folding scaffold found in other viral proteases [7, 8]. Domain 3 comprising of α-helices (residues 201–306 shown in red) is specific for CoVMpro. A loop of ~15 residues (185–200 shown in purple) links domains 2 and 3. A H41–C145 catalytic dyad is present at the interface between domains 1 and 2 which hydrolyze the peptide bond at specific cleavage sites in a series of steps. Imidazole group of H41 extracts the side-chain proton of C145 resulting into thiolate nucleophile which then attacks the amide bond of the substrate. As a result of proton extraction, N-terminal portion releases from H41 which in turn releases the C-terminal portion restoring the catalytic dyad. Therefore the inhibitors which perturb the communication between the H41 and C145 would block the catalytic cycle of the protein. N3 is a computationally designed peptido-mimetic irreversible Michael acceptor inhibitor which has the capability to inhibit Mpro of not only SARS-CoV2 but also other coronaviruses, including SARS-CoV and MERS-CoV [9]. It acts as a natural target of Mpro that forms a covalent bond with the catalytic site residue like any Michael acceptor inhibitor, blocking the catalytic cycle of the protein [10]. The inhibitor also showed antiviral activity against bronchitis virus in an animal model [11]. The 50% cytotoxicity concentration of the inhibitor is >133 μM [9].

Fig. 1

Structure of Mpro. Domain 1, orange; domain 2, green; domain3, red. Domains 2 and 3 are connected by an unstructured loop shown in purple. The catalytic site is magnified showing H41-C145 catalytic dyad

Developing a new drug for a disease is a time-taking process. While such processes of drug development are running, there is a need of rapidly finding a rescue for the interim time. Repurposing of pre-existing drugs for the still escalating disease which took more than 1.9 million lives as of January 13, 2021, around the globe is one of such faster and more efficient approaches [12-15]. Continuous endeavor for the development of effective treatment measures by researchers worldwide has suggested some potential drugs for repurposing both in silico and experimentally. In March 2020, the WHO commenced the “SOLIDARITY Trial” of four prevailing antiviral drugs for targeting different proteins of the virus on thousands of affected people in 10 countries to estimate the treatment effects of the drugs. The drugs chloroquine and hydroxychloroquine, ritonavir/lopinavir with and without interferon-beta, and remdesivir are presently under clinical trials of different phases [16]. Additionally, a number of antibodies, antivirals, and various therapeutic compounds have been reported [2, 12–15, 17–28]. However, till date, none of these drugs or compounds can be mentioned as part of the confirmed treatment protocol for the disease.

Computer-aided drug designing (CADD) offers methods to accelerate the search for possible lead compounds, and the virtual screening of compounds with respect to a receptor is one of such widely used methods [29-33]. The standard practice of such screening is to use a conformation of the receptor as a reference to find suitable ligands that fit in its binding pocket; the model of the conformation may come either from experiment or from computation or their combination. But in a real system, in principle, all possible conformations are present, and therefore a ligand can select the conformation of the receptor that it suits, even if that is not the most probable conformation of the uncomplexed/apo receptor, and then through a population shift, the equilibrium can gradually shift towards the complex state. Therefore, the use of one conformation of the receptor, even if it is experimentally determined, can miss the ligands that would be of higher binding affinity for another conformation of the binding pocket of the receptor. To overcome this limitation, a relatively better approach is to use multiple conformations of the receptors as the references for screening the ligands. In this work, starting from the experimentally determined structure, molecular dynamic (MD) simulation-based conformational sampling has been used to generate an ensemble of conformations, and from that pool, fourteen different conformations of the receptors along with the experimentally derived receptor were chosen for the screening using the “DrugBank” database containing 2115 FDA-approved compounds [34].

Methods

Molecular dynamics simulation

The uncomplexed (Apo) structure was modeled by removing the ligand N3 from the structure of the SARS-CoV-2 main protease (resolution 2.16 Å) having the PDB ID 6LU7 [9]. CHARMM36 force field [35] and parameter set were used. The system was solvated using a cubic water box of TIP3P [36] water model. K+ and Cl− ions were added to neutralize the overall charge of the system and also to maintain a 0.15 M ionic concentration. The energy minimization of the model and molecular dynamics (MD) simulation were performed using NAMD2.12 package [37], at 300 K, under NPT conditions, using a 2 fs time integration step, and the snapshots were saved at every 4 ps interval. Heating and equilibration were followed by production run. Short-range nonbonded interactions were truncated at 12 Å, whereas the long-range electrostatic interactions were computed using the particle mesh Ewald (PME) method [38]. The pressure and the temperature were controlled using Nosé-Hoover thermostat/barostat algorithm [39]. A 100 ns trajectory was run. Before continuing with the docking procedure, the trajectory was clustered skipping every alternate frame in UCSF Chimera based on backbone RMSD of the protein with default RMSD cutoff provided by the software [40]. Because of the relatively large population of the first 14 clusters, central frames of these clusters (named as frames 2–15) were selected. The apo crystal structure (named as frame 1) along with frames 2–15 were selected for ensemble docking. (Fig. 2, Fig. S1).

Fig. 2

(a) Differences in the side-chain orientation of some important residues present within the catalytic site shown with some representative snapshots chosen for screening (frame 3 and 9). The conformations are distinguished by the coloring schemes; among them, the pink color refers to 6LU7. (b) and (c) Interaction of N3 with Mpro in frames 3 and 9 respectively

Protein preparation

Apo structures, one extracted from 6LU7 and the others selected from the MD simulations, were prepared using the Protein Preparation Wizard of Schrödinger software suite 2019-1 release. Amino acid bond orders were adjusted following the addition of missing hydrogen atoms. Epik was used to calibrate the protonation as well as tautomeric states of Lys, Arg, Glu, Asp, and His to match pH 7(± 2). Finally the geometry was refined by restrained minimization with heavy atoms RMSD 0.3 Å.

Receptor grid generation

A 20 Å 3D grid box of the prepared receptors was defined keeping the important catalytic residues in the center. van der Waals radius scaling and partial charge cutoff were kept default.

Ligand preparation

The downloaded DrugBank database [34] and N3 ligand obtained from the complex were prepared using LigPrep without considering the protonation and ionization states. Also no stereoisomers were generated.

High-throughput virtual screening (HTVS) and docking

Using the prepared receptors and the ligands from the DrugBank, rigid HTVS was done in Glide [41, 42] version 8 with default settings. Glide uses an empirical scoring function which has been parameterized for separating binders from the nonbinders in screening. The best 25% of the HTV-screened ligands were used for SP (standard precision) docking. Further, the best 25% of SP-docked ligands were taken for XP (extra precision) docking [43]. In XP docking, a more sophisticated scoring function is included so as to wipe out the false-positive results which SP docking permits.

Binding free energy calculation

The free energy of ligand binding was calculated on the docked complexes using the Prime MM-GBSA method implemented in the Schrödinger software suite using default parameters. Prime MM-GBSA calculated the OPLS2005 force field–based energies of each system, using the VSGB2.0 implicit solvent model to account for the solvation energy [44].

Results and discussion

Selection of the MD frames for screening

The trajectory obtained after the simulation was first visualized in VMD, and no major deviation of the protein backbone from the starting conformation was observed which was also reflected in the root-mean-square deviation (RMSD) of Cα along the trajectory (Fig. S1). From the cluster of conformations, the central conformation of each of the top 14 clusters was selected for screening, in addition to the one obtained from the PDB structure itself. The Cα RMSD of each selected conformation, with respect to the PDB structure (PDB ID: 6LU7), is listed in Table 1. Although the overall backbone conformation of the protein remained stable, yet in the catalytic site significant variations in the amino acids side chain orientations were noticed (Fig. 2, S2).

Table 1

Cα RMSD (with respect to the 6LU7) of the conformations chosen for the screening. The ligand N3, as found from the crystallized complex (6LU7), was also docked to all other frames to get the reference score for all the frames

Frame index referred to in the screening	RMSD with respect to 6LU7 (Å)	XP g-score of N3 (kcal mol−1)	Binding free energy (kcal mol−1)
1( 6LU7)	0	−7.5	−82.4
2	1.324	−7.5	−78.9
3	0.998	−6.9	−75.7
4	1.122	−7.0	−93.2
5	1.338	−3.7	−74.8
6	0.857	−5.0	−74.5
7	1.143	−6.7	−83.1
8	0.097	−4.5	−70.6
9	1.324	−2.2	−70.0
10	1.452	−3.7	−48.2
11	1.326	−3.1	−81.4
12	1.507	−5.1	−66.1
13	1.848	−4.5	−78.4
14	2.516	−5.6	−72.8
15	2.185	−5.6	−71.4

Noting the interactions between Mpro and N3

N3 was previously reported to interact with T24–T26, H41, M49, F140, L141, G143-C145, H163–H164, M165, E166, P168, H172, and D187-Q192 [19]. In the crystal structure, N3 has been found to form a covalent bond with C145 [9]. Re-docking of N3 in all the 15 frames using the docking pipeline crossed all of the stages; among them, frame 1, which correspond to the receptor’s conformation seen in the crystal stature, showed the best XP g-score (−7.5 kcal mol−1), along with frame 2, as listed in Table 1. When frames 1 and 2 are further compared in terms of the binding free energy, frame 1 appeared as the best. Therefore, the positive control run successfully selected the best receptor’s conformation for itself, being fully in agreement with the experimentally determined conformation. Visually, it was observed that the binding pose of N3 was similar to that of crystal structure in multiple frames, and that serves as the benchmark of the method and the capacity to identify the best pair of receptor’s conformation and ligand (Fig. 3a). The g-score of the 6LU7 was used as a reference to find a ligand-drug with higher affinity. Ligand interaction diagram (Fig. 3b, c) of the re-docked complexes revealed the interaction of N3 with the catalytic dyad H41-C145 of the active site. These are the key residues for a proteolytic process where the thiol group of the C145 acts as nucleophile [45]. N3 was also found to form a hydrogen bond with N142, G143, E166, and Q189 and salt bridge with E166 in different frames (Fig. 3b, c). The salt bridge formation between N3 and Mpro is very important because E166 is known to have a crucial role in the dimerization of Mpro by associating the binding site with dimer interface, and that is necessary for the catalytic function of the protein [31].

Fig. 3

a) Superposition of N3 docked in the catalytic site in different frames. The conformations are distinguished by the coloring scheme. The crystalized conformationl of the receptor is shown in pink, and that of N3 is shown in green. (b) and (c) 2D diagrams of the interactions between Mpro and N3. The residues are coloured according to their nature of sidechain: hydrophobic (green) color, positively charged (blue), negatively charged (red), other polar groups (sky blue), and glycine (white)

Shortlisting the ligands

The top-scoring ligands obtained from XP docking were analyzed to find their type and strength of interactions. Such a blind screening listed 78 ligands from the entire database, which had scores better than N3 in at least one of the 15 frames of the receptor. Among them, there were a number of dyeing agents, food additives, supplements, and some nondrug molecules which were eliminated from the list; for obvious reasons, the drug molecules approved for topical application were also eliminated. This knowledge-based scrutiny yielded a shortlist of 49 compounds (Table S1, Fig. S3). Since the N3 itself showed a very similar or very close value of the dock scores in multiple frames, e.g., frames 1, 2, and 4, the short listed compounds which were present in at least 20% of the frames, i.e., 3 frames, were further identified (Table S2) which resulted a list 23 compounds. In this list, 9 more compounds were added, which otherwise got eliminated due to shortfall in the number of frames in which they appeared. The reason of adding those compounds in the list was that they were biochemically and/or functionally close to one or more of those 23 compounds. For example, doxorubicin and valrubicin (Table 1, Serial no. 13 and 30, respectively) were present in 3 frames each, but three other, rubicins amrubicin, daunorubicin, and idarubicin were left out because of their presence in less number of frames; they were added in Table 2 because of their biochemical and functional similarity with doxorubicin and valrubicin.

Table 2

Top-scoring ligands with their best XP g-score, number of frames in which they were found, and corresponding binding free energy

Serial no.	Name of drug	No. of frames where present	Highest g-score	Binding free energy in the highest-scoring frame (kcal mol−1)
1	Almasilate	6	−7.9	82.74
2	Amrubicin	2	−7.3	−74.27
3	Arbutamine	3	−8.9	−73.89
4	Cangrelor	11	−8.4	−47.83
5	Carbetocin	6	−8.0	−81.32
6	Cefixime	1	−7.8	−43.27
7	Cefpiramide	5	−8.1	−77.76
8	Ceftolozane	3	−8.9	−73.4
9	Daunorubicin	1	−6.8	−65.1
10	Desmopressin	11	−9.9	−59.36
11	Diosmin	9	−10.3	−55.58
12	Dobutamine	1	−8.2	−74.9
13	Doxorubicin	3	−7.7	−73.7
14	Droxidopa	5	−7.8	−44.7
15	Enviomycin	2	−7.8	−42.1
16	Felypressin	9	−8.9	−86.8
17	Idarubicin	1	−6.2	−56.9
18	Kanamycin	2	−7.9	−46.8
19	Lymecycline	4	−8.0	−73.31
20	Mitoxantrone	3	−9.3	−79.8
21	Novobiocin	2	−8.0	−74.3
22	Polydatin	5	−8.1	−66.0
23	Regadenoson	6	−7.7	−54.8
24	Ribavirin	4	−7.1	−29.9
25	Risedronate	6	−7.9	−21.34
26	Sofosbuvir	1	−7.3	−75.7
27	Terlipressin	3	−9.8	−84.2
28	Ritonavir	1	−7.6	−73.1
29	Tigecycline	3	−7.3	−67.2
30	Valrubicin	3	−8.1	−85.6
31	Zanamivir	1	−7.8	−46.0
32	Zoledronic acid	4	−7.7	−16.3

From this list, the ligands which were present in at least 3 frames were elaborated further in with the details of their interaction pattern with Mpro (Fig. S3). Interestingly for most of the compounds, the H-bond and salt bridge with E166 were common indicating that they might have the capability to disrupt the dimerization on Mpro which is important for performing the catalytic function [31]. A number of the screened drugs were found to interact with H41. The H-bonds and/or other interactions between H41 and different ligands could be argued to disrupt the interaction of the H41-C145 catalytic dyad rendering Mpro unable to perform its catalytic function.

Complexes from screening

The Mpro, which was derived from crystal structure 6LU7 by eliminating N3, identified ritonavir as a promising drug candidate although the g-score was lesser than a number of other hits. The score − 7.6 kcalmol−1 however was almost the same as N3; nevertheless, the identification of ritonavir provides a confidence to this computed prediction as the ligand is currently under discussion everywhere. The interaction pattern of ritonavir with residues of Mpro as depicted in Fig. 4 revealed this ligand also interacts with the catalytic site residues especially H41-C145 catalytic dyad forming π-π stacking interaction with H41 which might disrupt the communication between the dyad residues. In addition to ritonavir, a number of antivirals along with antibiotics, vasodilators, antineoplastic drugs like different rubicins, drugs used to treat circulatory diseases, drugs for bone diseases, etc. were found as top-scoring ligands in different frames (Table 2).

Fig. 4

Orientation and interactions of ritonavir in the catalytic site: (a) 2D ligand interaction diagram that shows stacking interaction between ligand and H41 (green), (b) 3D orientation of the ligand in the same conformation. The color scheme followed Fig. 1

Antivirals like ribavirin, sofosbuvir, zanamivir, penciclovir, saquinavir, ganciclovir, and cidofovir were found, some of them having known protease inhibitor properties. Most of the ligands were found in a single frame or in multiple frames with very low g-score (data not shown). All these drugs though have been found as potent Mpro inhibitors in other studies [46-51]; here, only ribavirin, which is presently under phase II clinical trial, was found in 4 frames, the g-scores being less than N3 and ritonavir (Table S3). The score of ribavirin was the highest in frame 6 (−7.179 kcal mol−1), and plausibly is due to the formation of H-bonds with residues R188 and T190 in this frame but not elsewhere (Fig. S5a). Further, in the lowest-scoring frame, i.e., in frame 13, the NH2 group of the ligand was not involved in the formation of H-bond with any residues of the catalytic site (Fig. S5b). Binding free energy for all the frames ranging between −29 and −39 kcal mol−1 was much lower compared to that of other ligands. The reason of a low g-score and binding free energy could be attributed to a smaller size of the ligand which leads to lesser accessibility to catalytic site residues.

Several antibiotics are presently being used to treat different viruses including SARS-CoV [52, 53]. Azithromycin is being used to treat COVID-19 but possibly to prevent secondary infections and may not be interfering with the components of the virus [53]. The antibiotics screened in this work, enviomycin, lymecycline, novobiocin, kanamycin, tigecycline, cefpiramide, ceftolozane, and cefixime, have been already known to be used as antivirals as well as for lung diseases. Lymecycline, which has been found as an inhibitor for Mpro in different studies, was identified in 4 frames. There was π-π interaction between lymecycline and H41only in the highest-scoring frame, i.e., frame 2 (Table S4, Fig. S6a). Also, H-bond and salt bridge are formed between E166 and NH3+ and NH2+ group of the ligand in this frame while in other frames salt bridge was formed with NH2+ group (Fig. S6b). The salt bridge between NH2+ group of the ligand and E166 might be considered as capable of stabilizing the ligand within the catalytic site of Mpro, but it can be stated that the involvement of both NH2+and NH3+in the formation of salt bridge as well as H-bond with E166 led to better stabilization as well as optimization of the ligand resulting the best g-score as well as binding free energy (Table S4). Cefpiramide, present in 5 frames, was found to have the best score in frame 7. In this frame, H-bonds were more in number. Additionally two salt bridges were formed with H41 and E47 (Fig. S6c). Tigecycline was found in 3 frames g-scores varying between −6.9 and −7.3 kcal mol−1. Salt bridge was present between the NH2+ group of the ligand and E166 in all frames. Only in the highest-scoring frame, there was a cation-π interaction with H163 (Fig. S6d). The highest score of ceftolozane though is better than N3 (−8.9 kcal mol−1); in other three frames, the scores were much less than N3, the lowest being −4.8 kcal mol−1. Only in the highest-scoring frame there was a salt bridge and π-π stacking interaction between H41 and the ligand (Fig. S6e). In lowest-scoring frame in addition to H41, salt bridge was formed between the NH3+ group of the ligand and E166 which might have over-stabilized the ligand resulting in a low g-score (Fig. S6f). Binding free energies for all the antibiotics were in agreement with the g-scores (Table S4). The rest of the antibiotics, though with g-scores better than N3, were found in less number of frames, and therefore they were not discussed here in detail.

Among the top-scoring ligands, rubicins, the anthracycline antineoplastic antibiotics were found in 6 frames. Previous studies have reported some of the resulted rubicins to have in vitro effect to the replication of other viruses though their protease-inhibiting capacity was not clear [54-56]. These drugs along with the detail of their interactions have been listed in Table S5. Most of the drugs had g-scores greater than or very similar to N3. Among the 3 frames in which doxorubicin was found, g-scores in frame 4 (−7.4 kcal mol−1) and 13 (−7.7 kcal mol−1) were mostly similar with the score of apo-N3 while in frame 8, it was less than N3 (−6.4 kcal mol−1) and it was seen in the ligand interaction diagram that in this frame, there were both cation-π and π-π interactions with H41 (Fig. S7b) while in the highest-scoring frame (frame 13), there was only π-π interaction with H41 (Fig. S7a). Also in the lowest-scoring frame, E166 was forming two H-bonds with the ligand in contrast to no H-bonds between E166 and the ligand in the highest-scoring frame. Valrubicin was found in frames 3, 7, and 13 among which scores in the first two frames were higher than N3 (−7.7 kcal mol−1and − 8.1 kcal mol−1, respectively) (Table S5). Only in the lowest-scoring frame (−5 kcal mol−1) π-π interactions with H172 were found (Fig. S7d). In both the above rubicins, binding free energies were in harmony with the respective g-scores (Table S5). Other rubicins, daunorubicin and idarubicin, were found in frame 4, but their g-scores were less than N3. Both of them were found to form salt bridges as well as H-bonds with E166 (data not shown).

One crucial reason for fatality in COVID-19 patients is venous thromboembolism [57-61]. The use of drugs that inhibit intravenous clotting like heparin and tissue plasminogen activator along with other drugs have already been disclosed as more effective [62-66]. In the present investigation, cardiovascular drugs like cangrelor and diosmin were found in more than 50% of the frames and imparted in more rigorous interactions in MD-generated receptor frames suggesting optimization of the complexes. Residues involved in the interactions with diosmin and cangrelor in different frames are furnished in Tables S6 and S6, respectively. Among the different plant flavonoids reported as SARS-CoV-2 inhibitor, diosmin has been reported in several studies [67-71]. In our multi-conformation-based screening study, diosmin, an anticoagulant, was found in 9 frames, the highest and lowest being frames 10 (−10.266 kcal mol−1) and 6 (−7.517 kcal mol−1), respectively (Table S6). In the ligand interaction diagram, it was seen that there were cation-π and π-π interactions among the ligand and H41 of Mpro in the lowest frame (Fig. 5c). Binding free energy of lowest-scoring frame was however more or less similar to that of the highest-scoring frame. The ligand orientation within the binding pocket varied in the highest and the lowest-scoring frames. L167 was forming H-bond with the ligand in the lowest-scoring frame. While g-scores of diosmin in most of the frames were more than or near the values of N3, scores of cangrelor, another cardiovascular drug, were mostly less than that of N3 (Table S7). In 2 out of 11 frames, the value was more than that of N3. H41 was found to form salt bridges with the molecule in all the frames except in the lowest one. The orientation of the ligand was different in the highest and lowest-scoring frames. In the highest-scoring frame, side chains E166 and K137 were oriented towards the outer side of the catalytic side, so they were not involved in the H-bonding (Fig. 6b). The PO4-containing group of the ligand was forming H-bonds with the catalytic site residues. Cl− of the ligand was also involved in the interaction in the highest-scoring frame. On the contrary, in the lowest-scoring frame, K137 was oriented towards the ligand and was forming salt bridge with the PO4 group (Fig. 6c, d). E166 was also involved in forming H-bond with the ligand. Regadenoson was screened in 6 out of 15 frames. Details of catalytic site residues involved in the interaction with this drug as depicted in Table S8 showed the g-scores were less than N3, the lowest and highest being −5.156 and -7.678 kcal mol−1, respectively. The number of H-bonds were more in the highest-scoring frame (frame 7), and the –OH groups of the ligand were involved in H-bond formation (Fig. 7a) while in the lowest-scoring frame (frame 11), the –OH groups were not forming H-bonds (Fig. 7c). Binding free energies of cangrelor and regadenoson were not consistent with the g-scores. It can be argued that the best pose of a ligand might not result in the lowest binding free energy as the type of interactions could alter the fitting of the ligand. Again, the poses with lowest binding free energy might be associated with over-stabilization of the ligands resulting in perturbation of the optimal functionality of the ligand.

Fig. 5

2D ligand interaction diagram and 3D orientation of the diosmin are shown in (a) and (b) respectively for the highest scoring frame. Cation-π (red) and π-π (green) interactions among the ligand and H41 of Mpro in lowest frame 3D orientation of the ligands, and functional groups involved in the interaction are different in the highest- (b) and lowest (b)-scoring frames. Color scheme followed Fig. 1

Fig. 6

2D ligand interaction diagram and 3D orientation of cangrelor are shown in (a) and (b) respectively for the highest scoring frame, which reveal the  H-bond between PO4 group of the ligand and the catalytic site residues. The same for the lowest scoring frame are shown in (c) and (d) which show the salt bridge between K137 and PO4 group.. Color scheme followed Fig. 1

Fig. 7

2D ligand interaction diagram and 3D orientation of regadenoson are shown in (a) and (b) respectively for the highest scoring frame. The same for the lowest scoring frame are shown in (c) and (d). The higher number of H-bonds and the involvement of –OH group in case of highest scoring frame compared to the lowest scoring frame are notable. Color scheme followed Fig. 1

Three vasodilators desmopressin, felypressin, and terlipressin, were screened as Mpro inhibitors among which desmopressin and felypressin were found in 11 and 10 frames, respectively. Though there was no experimental evidence of antiviral effects of the pressins, felypressin has been found as an inhibitor of the dengue virus in silico [72]. Details of interactions of desmopressin and felypressin are given in Table S9 and S10, respectively. Side-chain orientations of catalytic site residues of Mpro and also the ligand desmopressin in the highest (frame 4) and lowest (frame 11) scoring frames were different (Fig. 8). In frame 4, N119 and N142 were not involved in H-bonds while T190 and Q192 were forming H-bonds with the ligand (Fig. 8a, b). In this frame, T190 and Q192 were not oriented towards the ligand in frame 11 and thus were not involved in H-bonds formation. XP g-scores of felypressin were the highest and lowest in frames 4 and 13, respectively (Table S10). The number of H-bonds between ligand and Mpro was much higher in frame 4 especially with E166 (Fig. 9a, c). Also salt bridges were present between E166 of Mpro and one of the NH3+ of the ligand in frame 4. One of the NH2 instead of NH3+was forming H-bond with E166 in frame 13, and H41 was forming H-bond with the ligand (Fig. 9b, d). The ligand was nearer to side chain of C44 in frame 13 where there was an H-bond between C44 and ligand, but in frame 4, the ligand was away from C44, and also the orientation of C44 was different. The side chains of T190 and Q192 were oriented towards the ligand in frame 4 but were away from the ligand in frame 13. As a result, H-bonds between the residues were found in frame 4 but not in frame 13 (Fig. 9b, d). In the case of the highest-scoring frame of terlipressin (frame 1), like previous two pressins, one of the NH3+ group of the ligand was forming H-bond as well as salt bridge with E166 (Table S11, Fig. S9a). In the lowest-scoring frame, there was no salt bridge between the ligand and Mpro (Fig. S9b). Binding free energies were in correspondence with the g-scores.

Fig. 8

2D ligand interaction diagram and 3D orientation of desmopressin are shown in (a) and (b) respectively for the highest-scoring frame which show that T190 and Q192 are involved in H-bonds. The same shown for the lowest-scoring frame in (c) and (d) reveal that those H-bonds are missing. Color scheme followed Fig. 1

Fig. 9

2D ligand interaction diagram and the 3D orientation of felypressin are shown in (a) and (b) respectively for the highest-scoring frame which reveal the H-bond and salt bridge between NH3+ and E166. The same for the lowest-scoring frame are shown in (c) and (d) which find the H-bond formed between NH2+and E166. Color scheme followed Fig. 1

In addition to the above-discussed groups of drugs, a few drugs with different functions were found (Table 2). Some of them have been found as Mpro inhibitor in different studies [73-75]. Out of these drugs, almasilate, though found in 6 frames, is a buffering antacid; thereby, its chance of entering into the bloodstream is very less likely. Therefore, this drug has not been discussed in detail here. Also the binding free energies of the ligand were very low in different frames. Imami et al., with their transcriptome-based approach, have shown that carbetocin, the analogue of oxytocin, has an impact on the induction of immune response in COVID-19 [73]. Here, carbetocin was found as an Mpro inhibitor in 6 frames (Table S12) among which in the highest-scoring frame, the number of protein-ligand H-bonds was much more than that in other frames which might have resulted in the lowest binding free energy. Risedronate and zoledronic acid, drugs for the treatment of osteoporosis and other bone diseases, were found in 4 and 6 frames, respectively. Both of the ligands possess O− groups which were found to be involved in the formation of salt bridge with H41 (Table S12). In the case of zoledronic acid, additionally cation-π and π-π stacking interactions with H41 were seen. The binding free energy of these ligands was very low compared to that of others ranging from ~ −2 to −21 kcalmol−1. Polydatin, a picecid, used in various diseases like anti-arrhythmia drug and potassium channel modulator, was found as a probable Mpro inhibitor in previous studies [74, 75]. Here the ligand was found to be buried in the catalytic pocket of Mpro in the highest-scoring ligand while in the lowest-scoring one, a smaller number of residues were accessible (Table S12). Arbutamine and dobutamine, the drugs to treat cardiogenic shock, were found in frames 1 and 3, respectively. Low scoring frames of arbutamine showed covalent interactions with E166 which were not found in the highest-scoring frame. Rather in the highest-scoring frame, there was a π-π stacking interaction with H41. The binding free energy was highest in this frame. H41 of Mpro was found to form a π-π stacking and/or cation-π interaction, and E166 was involved in salt bridge formation with droxidopa, a norepinephrine prodrug, in all frames except the highest-scoring one (frame 6). Mitoxantrone, an antineoplastic drug, was obtained in 3 frames among which a number of protein-ligand H-bonds were least in the highest-scoring frame (Table S12). The binding free energies of droxidopa and mitoxantrone were in harmony with the g-scores. It can be stated that more protein-ligand interaction might sometimes result in over-stabilization of the ligands which might alter the optimal function of the ligand. Some of the drugs discussed in this paragraph were detected as an inhibitor of Mpro in previous studies, most of them being in silico. These drugs could be tested further experimentally. Modification of functional groups of the rest of the drugs might result into new compounds with the capability to inhibit Mpro.

The dock score of N3 was taken as reference, which is experimentally known to bind to Mpro and used it to shortlist the hits with scores better than N3. This was done with the expectation that it would be able to identify the ligand at least having similar interaction pattern with desired biochemical activity and then would carry some additional features that might increase the probability of success rate of the prediction. The present virtual screening and MM-GBSA-based approaches have revealed a number of FDA-approved drugs which have Mpro affinity above the cutoff of N3 scores. Almost all the top-scoring drugs were forming hydrogen bonds with the catalytic site residues H41, M49, F140, L141, G143-C145, H163–164, M165, E166, P168, H172, and D187-N192, the same contact points seen in the N3 complex. Though not all, some of the identified drugs have already proven antiviral effects.

Variation in the g-scores and binding free energy for the same ligand in different frames were due to differences in the interactions and packing of the ligand inside the catalytic site. There have been some exceptions in the trends of the g-scores and the binding free energies, which otherwise reflected the correlation in their trends, i.e., the good g-score did not correspond to promising binding affinity in some cases which indicated the need of some advanced method of refinement. But in this work, shortlisting of ligands by either of the scoring methods was the goal, and as that was achieved, the strict match between g-score and binding free energy was not aspired. Again the consonance of binding free energy with the g-score in most of the ligands has confirmed that the particular pattern of interaction in the best scoring ligands resulted in the lowest binding free energies. Also, it was evident that the usage of multiple conformations for screening resulted in some molecules which might have been rejected if a single conformation would have been used because of relatively lower g-score and binding free energy.

Conclusions

A list of FDA-approved drugs have been obtained through virtual screening with the promise to inhibit Mpro of COVID-19. The N3 molecule is already known to bind to Mpro, which is reported as its co-crystallized ligand. The binding score of N3 has been used as a reference, and the ligands with scores better than this are plausibly expected to bind to Mpro with higher affinity. The scores of the docking and also the computed binding free energies are best interpreted in a relative scale, instead of debating on the absolute values. The use of ensemble representation of the receptor’s conformation has implicitly covered the possibilities of “induced-fit” and/or “population shift” mechanisms that the ligands might adapt to fit in the pocket, and therefore this method has been powerful to identify a wide variety of ligands, and many of which would have been otherwise missed out. It has not only identified antiviral drugs like ribavirin, ritonavir, etc., but it has also identified the drugs used in various other diseases, e.g., amrubicin, cangrelor, desmopressin, diosmin, etc. as the potential candidates to inhibit Mpro. The screened drugs might have a cumulative inhibitory effect on the viral multiplication and pathogenesis. Functional group modification of the drugs might lead to better binding affinity with Mpro. Top-scoring ligands with correlating binding free energies could be considered for persuasion of further tests and designing inhibitors against the viral protein.

(PDF 989 kb)