3CLpro and PLpro affinity, a docking study to fight COVID19 based on 900 compounds from PubChem and literature. Are there new drugs to be found?

The spread of a novel coronavirus SARS-CoV-2 and a resulting COVID19 disease in late 2019 has transformed into a worldwide pandemic and has effectively brought the world to a halt. Proteases 3CLpro and PLpro, responsible for proteolysis of new virions, represent vital inhibition targets for the COVID19 treatment. Herein, we report an in silico docking study of more than 860 COVID19-related compounds from the PubChem database. Molecular dynamic simulations were carried out to validate the conformation stability of compound-ligand complexes with best docking scores. The MM-PBSA approach was employed to calculate binding free energies. The comparison with ca. 50 previously reported potential SARS-CoV-2's proteases inhibitors show a number of new compounds with excellent binding affinities. Anti-inflammatory drugs Montelukast, Ebastine and Solumedrol, the anti-migraine drug Vazegepant or the anti-MRSA pro-drug TAK-599, among many others, all show remarkable affinities to 3CLpro and with known side effects present candidates for immediate clinical trials. This study reports thorough docking scores summary of COVID19-related compounds found in the PubChem database and illustrates the asset of computational screening methods in search for possible drug-like candidates. Several yet-untested compounds show affinities on par with reported inhibitors and warrant further attention. Furthermore, the submitted work provides readers with ADME data, ZINC and PubChem IDs, as well as docking scores of all studied compounds for further comparisons.

Introduction

In the late 2019, the world was taken by storm by a novel type of β coronavirus dubbed SARS-CoV-2 or 2019-nCoV originating in Wuhan, China [1,2]. The disease caused by SARS-CoV-2 was named COVID19 and like diseases caused by other coronaviruses outbreaks, such as SARS-CoV [3] (Severe Acute Respiratory Syndrome Coronavirus) in 2003 and MERS [4] (Middle East Respiratory Syndrome Coronavirus) in 2012, it presents itself with respiratory symptoms including fever, dry cough, fatigue and loss of taste or smell. In more severe cases COVID19 leads to pneumonia, dyspnea and eventually death [5,6].

Key moments for stopping the spread of a virus in the human body, besides vaccination, include: A. virus entrance to the cell; B. inhibition of polyproteins proteolysis, whose products are subsequently employed in new virions production; C. replication of the RNA genome. Thus, potential targets for SARS-CoV-2 treatment include: A. the SARS-CoV-2 receptor-binding spike protein; B. proteases 3CLpro (also denoted as Mpro)/PLpro (3C-like protease/papain-like protease); C. the RNA polymerase (RdRP) [7,8]. SARS-CoV-2 spike protein binds angiotensin-converting enzyme 2 (ACE2) making it an ideal target for vaccine development, preventing the virus from entering the cell. Proteases 3CLpro/PLpro and RNA polymerase display high genome sequence similarities with SARS-CoV (82–96%) and present inhibition options A. and B. for terminating the virus replication in cells [9]. Protease 3CLpro appears to be the most promising inhibition target of SARS-CoV-2 polyproteins proteolysis (B). 3CLpro monomer consists of three domains (domain I – residues 8–101, domain II – residues 102–184, domain III – residues 201–303) and a long loop (residues 185–200) connecting domains II and III.

Its active site forms a cavity located in the cleft between domains I and II with catalytic dyad (Cys145 and His41) that is large enough to bind the active substance forming a drug-protein complex [10,11], see Fig. 1 A. The interior of the cavity is composed of several hydrophobic amino acids (Leu27, Phe140, Gly143, Met165 Pro168), forming hydrophobic environment and polar amino acids (Thr26, Asn142, His163, Gln189) capable of stabilizing the drug-protein complex [10,11]. 3PLpro consists of two distinct domains (ubiquitin-like domain – residues 1–60, thumb-palm-finger domain – residues 61–315) with its active site containing a catalytic triad (Cys111, His272, Asp286) and an important zinc ion coordinated by four cysteine residues [12], see Fig. 1B. The surrounding of a smaller catalytic site is comprised of polar amino acids (Asn109, Asn110, Thr265, Tyr273) and hydrophobic amino acids (Gly266, Gly271, Gly287, Ala288).

Fig. 1

Three-dimensional ribbon structures of 3CLpro (1A) and PLpro (1B) proteases, with catalytic diad/triad indicated by sticks and balls representations of participating molecules.

In the herein presented in silico work, we have performed molecular docking experiments of potential SARS-CoV-2 inhibitors accessed from external compounds database against two 3CLpro and PLpro structures. Docking protocols were followed by molecular dynamics simulations of most promising compound-protein complexes to ascertain their conformational stability, and by the evaluation of free energy of binding. Several in silico research papers have been published over the previous year focusing either on docking studies [11,[13], [14], [15], [16], [17], [18], [19], [20], [21], [22]] of selected compounds or screening a large number of compounds from external databases utilizing machine learning approaches [23] or virtual screening [24]. Our obtained results were compared with previously published papers [11,[13], [14], [15], [16], [17], [18], [19], [20], [21], [22]] and redocking of several compounds with high affinity towards SARS-Cov-2 proteases from these works was performed to allow for a direct comparison to the presented results.

Materials and methods

Preparation

The 3-dimensional structures of 866 compounds tagged by the keyword “SARS-CoV-2″ [25] were downloaded from the PubChem [26,27] [October 23, 2020] database in the .sdf format and converted to the .mol2 format using the OpenBabel compiler 2.3.2 [28]. Furthermore, 3D structures of 33 best scoring compounds from four different publications [11,16,17,21] were downloaded from the PubChem database, ten best scoring compounds from another paper [13] were downloaded from the ZINC database [29,30] and three geometries of compounds not found in either database were optimized at the DFT level of theory (B3LYP/6–311G*) [31], [32], [33], [34], [35] using the Gaussian09 software [36]. Where possible, compounds were described by their CID number, DB identifier, InChI key and ZINC code as can be found in the ESI files. Polar hydrogens and Gasteiger charges were added to all structures using the AutoDock utility scripts [37].

The structures of SARS-CoV-2 3CLpro and PLpro at 100 K (PDB IDs: 6LU7 [38] and 6WZU [39]) as well as their crystal structures determined at the room temperature (PDB IDs: 6WQF [10] and 7CMD [40]) were downloaded from the RCSB protein data bank [41]. Swiss-PDBViewer [42] was used to fix the missing atoms within the proteins and to remove present ligands (6LU7 and 7CMD). Additional protein strands were removed from the selected crystal structures and the 6LU7 and 6WZU/7CMD structures were stripped of water molecules, while one crystal water was retained within the 6WQF structure. This molecule is situated within the pocket of the protease active site cavity, near residues His41 and Asp187, and it has been previously reported that this water molecule takes part in several interactions leading to charge stabilization on neighbouring residues [10].

Conformational differences between pairs of protein structures were evaluated by template modeling (TM) score ranging from 0 to 1, with 1 indicating a perfect match and 0 indicating no match at all between compared structures [43]. The selected 3CLpro structures have the TM score 0.98, with minor local structural differences at the C-terminus. The TM score further increases to 1 if calculated for the initial 300 amino acids only. The compared PLpro structures have the TM score 0.96, with minor differences throughout the whole 3D structures. Furthermore, 7CMD does not have a completely solved crystal structure, as there are amino acids missing between Val220 and Thr231. This should affect neither TM score, as it is calculated only from pairs of amino acids that are present in both structures nor the docking results, as the sequence is not located in the proximity of the active site cavity.

Docking

Semi-flexible docking simulations were performed using the software Autdock4.2.6 [37,44]. A grid box of 90×90×90 Å [with a resolution of 0.275 Å and centered at x, y, z = (-20 Å, -5 Å, 15 Å)] and a grid box of 60×60×60 Å [with same resolution and centered at x, y, z = (-8 Å, 82 Å, 37 Å)/(-32 Å, -24 Å, -28 Å)] were used to calculate potential maps within the protein binding cavity of 6LU7/6WQF and 6WZU/7CMD protein units, respectively. In the first approximation, total of 50 Lamarckian genetic algorithm docking runs were performed for each tested molecule, with the initial population containing 300 individuals. The maximum number of energy evaluations was set to 30,000,000 and the maximum number of generations was set to 27,000. Rates of gene mutations and crossovers were set to 0.02 and 0.8, respectively, with one individual to survive to the next generation. The probability of Solis & Wets local search was set to 0.06 with a maximum number of iterations of 300 and a minimum step size of 1% of the grid box. Compounds with docking scores below -11.0 kcal/mol (total of 34 structures) for 6QWF, below -12.0 kcal/mol for 6WQF (total of 63 structures) and below -9.0 kcal/mol for 6WZU/7CMD (total of 47/37 structures) were selected into a second round of docking calculations with the number of the genetic algorithm runs increased to 100 and the number of energy evaluations increased to 50,000,000. All resulting docked poses were clustered with 2.0 Å tolerance and analysed with the AutoDockTools utility [37,45] with one compound (CID: 73774610) omitted from further data evaluation due to unforeseeable errors.

Molecular dynamics (MD)

In addition, five structures with the lowest docking scores for 6WQF and 6WZU protein units were selected for MD simulations in GROMACS2018.7 [46], [47], [48], [49] using the CHARMM36-july forcefield and the molecular parameters and their topologies were generated with the CgenFF utility [50], [51], [52], [53], [54].

The proteins were solvated in the rhombic dodecahedron with edges of the box distant 1.0 nm with respect to the nearest protein residue using the four-point TIP4P rigid water model [55] and Na+/Cl− ions to ensure electro neutrality. The steepest descent algorithm was employed for energy minimization purposes, with the studied systems reaching their energy minimum within 1500 steps. Subsequently, all systems were equilibrated to a 300 K temperature (using the modified Berendsen thermostat [56]) and a 1 bar pressure (the Parrinello-Rahman coupling [57]) for a total of 200 ps. Finally, a 10,000 ps molecular dynamics simulation was run using the leap-frog integrator with a step of 0.002 ps. Resulting trajectories were analysed using the GROMACS software package calculating parameters such as the short-ranged compound-protein Coulombic interaction energy, the short-ranged compound-protein Lennard-Jones energy, the RMSD (Root-mean-square Deviation) of compound's heavy atoms with respect to the initial structure of the residues’ backbone, distances/angles between compound atoms and protein atoms within the cavity, etc. The formation and propagation of hydrogen bonds were analysed using the VMD 1.9.3 software [58] with a 3.5 Å distance cut-off and a 30° angle cut-off.

Free energies of binding were calculated using the MM-PBSA (Molecular Mechanics – Poisson-Boltzmann Surface Area) approach as implemented in the g_mmpbsa GROMACS package [59,60]. Snapshots at an interval of 10 ps were extracted from MD trajectories and binding free energies of the compound-protein complex were expressed as a difference between ensemble averaged free energies of complexes and their constituents.

The polar solvation calculation was carried out with a 0.1 M salt concentration and a solute dielectric constant of 2. The SASA (Solvent Accessible Surface Area) model was used in calculating non-polar contributions to free energies of binding.

ADME

Absorption, distribution, metabolism, excretion (ADME) parameters as well as other physicochemical descriptors and druglike properties of studied compounds were computed from a list of SMILES codes using the SwissADME website [61] and can be found in the supporting excel file PubChem_SwissADME.xlsx along with the other identifiers. These data were used to check whether the compounds satisfy criteria of empirical drug likeness rules, such as the Lipinski's et al. rule of five [62] and the Ghose et al. filter [63].

Results

Molecular docking analyses are nowadays an important part of drug design and are often used to examine drug-protein interactions to allow for a better understanding of the target active site or the importance of functional group substitutions within a drug structure. Final docked poses can be evaluated by the number of drug-protein interactions (H-bonds, π-π, etc.), as well as by the binding energy (score) which allows for a comparison of the affinity of a set of different compounds to the targeted receptor. Compounds with lowest (most negative) score and a reasonable hydrogen bond pattern can be chosen as perspective drug targets for further studies and derivatizations. Compounds downloaded from the PubChem database achieved in the initial docking runs against 6LU7 (6WQF) docking scores ranging from -2.80 (-3.35) kcal/mol to -13.88 (-14.45) kcal/mol, with a median score of -7.45 (-8.09) kcal/mol and an average score of -7.65 (-8.48) kcal/mol. Most negative docking scores were leaning heavily towards molecules with a higher number of atoms (molecular weight), which can be explained by the size of 3CLpro cavity. An average shift of -0.83 kcal/mol towards lower docking scores was observed, see Fig. 2 , when comparing the 6WQF docking scores against 6LU7 ones, with one outlier (CID: 3793) and a variance 0.45 kcal/mol of the difference between 6WQF and 6LU7 docking scores. The strong correlation of the docking scores is consistent with structural similarity as discussed in the previous section.

Fig. 2

The correlation of compounds’ docking scores for 6WQF 3CLpro to their respective 6LU7 3CLpro scores.

Docking in silico experiments with PLpro structure 6WZU (7CMD) yielded docking scores ranging from -2.77 (-0.72) kcal/mol to -10.51 (-11.21) kcal/mol, with a median score of -6.42 (-5.91) kcal/mol and an average score -6.61 (-5.98) kcal/mol with significantly lower bias towards compounds with a higher molecular weight. The average shift of -0.63 kcal/mol towards lower docking scores against PLpro structure 6WZU was observed when comparing to the respective docking scores of 7CMD, see Fig. 3 . The differences in selected PLpro crystal structures lead to a higher variance (1.91 kcal/mol) of the difference between 6WZU and 7CMD docking scores compared to the 6WQF and 6LU7 one (compare Figs. 3 and 2, respectively). It is evident that the origin of the targeted crystal structure affects docking score in the semi-flexible docking protocols. Due to a smaller cavity size, the differences in crystallization methods (crystallization temperature and/or liganded vs. un-liganded structure) affect the active site of PLpro, and hence the docking scores, more significantly than for 3CLpro.

Fig. 3

The correlation of compounds’ docking scores for 7CMD PLpro to their respective 6WZU PLpro scores.

See section “Comparison with other studies” for graphical representations of all obtained docking results. The docking studies against all four protein targets have not indicated any clear pattern between docking scores and the number of formed hydrogen bonds, as well as any of the most common ADME descriptors, with most complexes stabilized by six and less hydrogen bonds, see Fig. 4 . Redocking of top structures lead to a small decrease in docking scores indicating that the initially chosen docking parameters were sufficient for this study. Table 1 contains structures, CID identifiers, trivial names and docking scores of top five compounds studied along with the lowest docking score from the second round of redockings for all four protein targets. Extended version of this table containing top ten compounds for each protein can be found in Table S1.

Fig. 4

Frequency of predicted H-bonds formed between docked compounds and the respective protein targets.

Table 1

CID PubChem identifiers, trivial names (in brackets), docking scores in kcal/mol, recognized pharmacological functions and molecular structures [74] of five compounds with highest binding affinities towards 3CLpro/PLpro structures. Molecular structures are downloaded from the PubChem database [26,27]. For an extended version of this table see Table S1.

№	6WQF (3CLpro)	6LU7 (3CLpro)	6WZU (PLpro)	7CMD (PLpro)
1	5281040 (Montelukast)	3191 (Ebastine)	16923 (Solumedrol)	121893 (Protoporphyrinogen IX)
	-14.75	-14.11	-10.89	-11.58
	Anti-inflammatory and bronchodilating activity	Anti-inflammatory activity	Anti-inflammatory and immunosuppressive activity	–
	Image, table 1	Image, table 1	Image, table 1	Image, table 1
2	53472683 (Vazegepant/zavegepant)	53472683 (Vazegepant/zavegepant)	73425380 (TAK-599)	9548902 (Taurocheno-deoxycholate(1-))
	-14.57	-13.80	-10.69	-11.24
	Anti-migraine activity	Anti-migraine activity	Anti-MRSA activity	Human metabolite
	Image, table 1	Image, table 1	Image, table 1	Image, table 1
3	6918155 (Ciclesonide)	154573806 (GRL-024–20)	25245769 (Biliverdine(2-))	60947 (Tirofiban)
	-14.39	-13.58	-10.67	-10.97
	Anti-inflammatory and anitiviral activity	- [a]	Human metabolite	Anti-coagulant activity
	Image, table 1	Image, table 1	Image, table 1	Image, table 1
4	5459840 (20-Hydroxyecdysone)	5281040 (Montelukast)	135483998 (5-Methyltetrahydrofolate)	5362119 (Lisinopril)
	-14.07	-12.83	-10.66	-10.90
	Protective role in the cardiovascular system	Anti-inflammatory and bronchodilating activity	Anti-neoplastic and antidepressant activity	ACE inhibitor with anti-hypertensive activity
	Image, table 1	Image, table 1	Image, table 1	Image, table 1
5	154573806 (GRL-024–20)	46700782 (Razuprotafib)	122146 (-[b])	72734520 (Avitinib)
	-14.05	-12.83	-10.61	-10.77
	- [a]	Potential vasculature stabilizing activity	- [a]	Tyrosine kinase inhibitor
	Image, table 1	Image, table 1	Image, table 1	Image, table 1

(a) No pharmacological functions have been reported.

(b) Compound does not have a trivial name.

6WQF

The obtained screening results show that compounds such as the potential inhibitor of SARS-CoV-2 main proteinase 3CLpro Peptidomimetic aldehyde 11a [9] (-13.48 kcal/mol), the cathepsin K inhibitor Relacatib [64] (-13.49 kcal/mol), the pair of vitamin D analogues Calcitriol and Calcifediol (-13.63 and -13.78 kcal/mol, respectively), as well as the novel drug designed to treat acute respiratory distress syndrome Razuprotafib [65] (-13.86 kcal/mol) exhibit high affinities to 3CLpro (6WQF). Thus, indicating their potential in the treatment of SARS-CoV-2. The actual top five structures with the lowest docking scores against 6WQF obtained in this study include the compound GRL-024–20, with docking score -14.05 kcal/mol, whose crystal structure in the complex with 3CLpro protease had already been determined [66]. Another compound with a low docking score is 20-hydroxyecdysone, docking score -14.07 kcal/mol, which is also the only compound that is currently not part of any COVID19 treatment trials. Anti-asthmatic drug Ciclesonide [67] achieved docking score -14.57 kcal/mol with predicted H-bonds to Thr26 and Gln189 amino acids. Potential anti-migraine drug (phase three of clinical trials) Vazegepant/zevagepant [68], -14.57 kcal/mol, also emerged as a possible 3CLpro inhibitor with predicted H-bonds to Thr26, Asn142 and Gly143 amino acid residues. The lowest docking score against 6WQF, -14.75 kcal/mol, was observed for the anti-asthmatic drug Montelukast [69], see Fig. 5 A and B. Its predicted H-bonds pattern includes bonds between His163 and the carbonyl group of Montelukast as well as bonds between both atoms of its’ hydroxyl groups and Thr190/Arg188 residues in 6WQF cavity. Other notable frequently predicted H-bonds to 6WQF structure include bonds to Thr26, Gly143, Glu166 and Gln189 amino acids.

Fig. 5

Putative bindings sites of Montelukast, Ebastine, Solumedrol and Protoporphyrinogen IX in 6WQF, 6LU7, 6WZU and 7CMD protein cavities with the compounds in purple color and proteins in green, orange, blue and yellow, respectively. (A) Interacting amino acid residues of 6WQF with Montelukast. (B) Binding site of Montelukast in 6WQF 3CLpro. (C) Interacting amino acid residues of 6LU7 with Ebastine. (D) Binding site of Ebastine in 6LU7 3CLpro. (E) Interacting amino acid residues of 6WZU with Solumedrol. (F) Binding site of Solumedrol in 6WZU PLpro. (G) Interacting amino acid residues of 7CMD with Protoporphyrinogen. (H) Binding site of Protoporphyrinogen in 7CMD PLpro (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

6LU7

Screening results of the set of studied (PubChem chosen) molecules against the second 3CLpro structure 6LU7 differ from the 6WQF ones in relative positions of high affinity molecules in the top scores list. The lower docking scores to 6LU7 were observed in aforementioned Peptidomimetic aldehyde 11a [9] (-12.50 kcal/mol), anticoagulant Argatroban [70] (-12.76 kcal/mol), broad-spectrum anti-viral protease inhibitor GC-376 (-12.79 kcal/mol) currently under study for its potential for inhibiting SARS-CoV-2 main protease 3CLpro [71], and two compounds whose biological effect has not been described yet, CIDs: 122,146 and 71,316,139 (-12.57 and -12.66 kcal/mol, respectively). One of these compounds (CID: 122146) has also shown high affinity towards PLpro (see below). Five compounds with the lowest docking scores obtained in this study include previously mentioned Razuprotafib [65], with docking score -12.83 kcal/mol, Montelukast [69], score -12.83 kcal/mol, GRL-024–20 [66], score -13.58 kcal/mol, and Vazegepant (zevagepant) [68], with score -13.80 kcal/mol. It should be noted that these molecules exhibit different poses and H-bonds patterns than in the complex with 6WQF, see ESI Table 1A and B. The highest apparent affinity towards 6LU7, -14.11 kcal/mol, was observed for H1 antihistamine Ebastine [72] as already found in Vatansever et al. [73], see Fig. 5C and D. The predicted binding pattern of Ebastine consists of two hydrogen bonds between its oxygens and hydrogens of His163 and Glu166 amino acids. In a similar way to previous 3CLpro structure 6WQF, H-bonds with highest incidence include bonds from Thr26, Glu166 and Gln189 amino acid residues.

6WZU

Vazegepant/zevagepant [68] also emerged as one of compounds with highest binding affinity towards PLpro (-10.21 kcal/mol), along Taurochenodeoxycholate (1-) anion (-10.37 kcal/mol), taurochloric acid (-10.40 kcal/mol), immunosuppressant Hydrocortisone hemisuccinate (-10.41 kcal/mol) and folic acid (-10.58 kcal/mol). Top five best scoring compounds include human metabolite Biliverdine(2-), score -10.67 kcal/mol, with predicted H-bonds to Lys105, Trp106, Asp108 and Ala288 amino acids and the novel prodrug TAK-599 designed to treat MRSA (Methicillin-resistant Staphylococcus aureus) infection [75], score -10.69 kcal/mol, with predicted H-bonds to His89, Lys105 and Asp108 amino acid residues. The lowest docking score was observed for the anti-inflammatory immunosuppressant Solumedrol, see Fig. 5E and F, with predicted H-bonds to Lys105, Trp106, Asp286 and Ala288 amino acids. Excellent docking score, -10.61 kcal/mol, was also observed for compound with CID: 122146. This molecule achieved very good scores in both docking rounds against 6WZU and is the prime example of the importance of docking studies, with a good conformational stability within the protein cavity and a stable H-bonds pattern (see the next section) warranting further research. Overall, the predicted H-bond pattern of studied compounds in the complex with PLpro forms a smaller set of interactions found with a higher frequency, most notably H-bonds between Lys105, Trp106, Asp108, Ala288 amino acid residues and compound's carbonyl groups.

7CMD

The compounds with the lowest docking scores against PLpro 7CMD crystal structure significantly differ from results obtained for 6WZU structure. The only compound that is found among ten best scoring compounds against both PLpro targets is the prodrug TAK-559 [75] (-10.56 kcal/mol against 7CMD). Other low scoring compounds against 7CMD include yet not-studied compound with CID: 9810132 (-10.55 kcal/mol), angiotensin II antagonist Valsartan (-10.67 kcal/mol), human metabolite Geranylgeranyl diphosphate (-10.71 kcal/mol) and tyrosine kinase inhibitor Mastinib (-10.75 kcal/mol). Five best scoring compounds include another tyrosine kinase inhibitor Avitinib (-10.77 kcal/mol) that has been investigated for use in the treatment of non-small cell lung cancer (NSCLC) and its ability to supress cytokine storm associated with COVID19 made it part of clinical trials for the treatment of COVID19 [76]. Angiotensive enzyme (ACE) inhibitor Lisinopril (-10.90 kcal/mol) has also been already linked with the treatment of COVID19. It has been found that patients with prescribed ACE inhibitors or angiotensin receptor blockers (ARB) are at reduced risk of COVID19 while not increasing risk of ICU care [77]. Anticoagulant Tirofiban (trade name AGGRASTAT) (-10.97 kcal/mol) has also shown promise in reducing thrombotic effects observed in patients suffering from COVID19 and is part of further clinical trials [78]. Its predicted H-bonds include bonds with two amino acids from PLpro catalytic triad (His272 and Asp286), as well as bonds with Cys270 and Ala288. The second best scoring compound against 7CMD PLpro structure was the human metabolite Taurochenodeoxycholate(1-) (-11.24) formed by deprotonation of taurochenodeoxycholic acid. Its predicted H-bond pattern includes bonds with Lys105, Trp106, Asn267 and Cys270. However, its utility in treatment of COVID19 is unlikely (being human metabolite). The compound with highest apparent binding affinity towards 7CMD structure PLpro was Protoporphyrinogen IX (-11.58 kcal/mol), see Fig. 5G and 5H. This compound is a direct precursor of protoporphyrin IX, which plays a critical role in living organisms during production of hemoglobin or chlorophyll. Protoporphyrin has already emerged as a potential antiviral drug for treatment of COVID19 by method interfering with the interaction of ACE2 and the receptor-binding domain of spike protein of SARS-CoV-2 [79]. Its predicted H-bonds include bonds with Lys105 and Trp106. Most frequently predicted H-bonds in compound-protein (7CMD) complexes include bonds with Lys105, Trp106, Asp286 and Ala288.

Docked poses of abovementioned compound-protein complexes, together with the description of the predicted H-bonds can be found in ESI Table S2A–D. The docking results together with evaluation of drug likeness criteria can be found in ESI excel files Docking_6WQF/6LU7/6WZU/7CMD.xlsx.

Molecular dynamics – docking pose verification

To investigate the conformational stability and to overcome the rigid protein picture, molecular dynamic simulations in water solvent at physiological conditions were carried out on five complexes for the 6WQF 3CLpro structure (CIDs: 5,281,040, 53,472,683, 6,918,155, 5,459,840, 154,573,806) and for the 6WZU PLpro structure (CIDs: 16,923, 73,425,380, 25,245,769, 135,483,998, 122,146) for the time of 10 ns.

The conformational stability of formed compound-protein complexes was validated by calculating the relative RMSD value of compounds’ heavy atoms with respect to the initial structure of the residues’ backbone, see Fig. 6 A and B. Time evolution plot of RMSD for complexes with 3CLpro (6WQF) shows that four structures achieved equilibrium at the start of simulation with averaged RMSD values around 0.3 Å, with Montelukast (CID: 5281040) and Ciclesonide (CID: 6918155) manifesting small conformational changes at 3.2 and 4.5 ns, respectively, suggesting small drifts in the trajectory during the simulation. These values indicate that all four structures are well positioned within the cavity's active site and form stable complexes.

Fig. 6

MD time plot evolutions of RMSD values of compounds’ heavy atoms with respect to the initial structure of the residues’ backbone of 6WQF (A) and 6WZU (B).

RSMD values of three compounds’ heavy atoms in complexes with PLpro (6WZU) show a higher motion patterns during the simulation period. Based on the visualizations of dynamics trajectories and on the time evolution plot of hydrogen bonds, this increase in the RMSD value is caused by the movement of parts of molecules that do not participate in the formation of H-bonds while rest of the molecule remains firmly attached to the protein. On the other hand, 20-hydroxyecdysone (3CLpro) and Solumedrol (PLpro) have shown no conformational stability in the cavity and a gradual increase in the RMSD value, linked to a decrease in the compound-protein interaction energy, eventually being dislodged from the target proteins. Compared to other high scoring compounds, 20-hydroxyecdysone contains only a single carbonyl group and oxygens participating in formations of H-bonds identified by docking protocols belonged to hydroxyl moiety making it highly hydrophilic. Solumedrol, on the other hand, contains an ideal number of both, H-bond donor and acceptor atoms, and the reason for its compound-protein complex instability is unclear and requires further data collection.

A second parameter that was used to monitor the conformational stability of these complexes, as well as to examine the compound-protein interactions was the evolution of the number of hydrogen bonds formed between the target proteins and the studied compounds during the MD simulation. The MD simulations have shown that the 3CLpro complexes form up to five hydrogen bonds, but these do not appear to be long-lived and the actual number of stable H-bonds is between one and three. In comparison, complexes with PLpro exhibit formation of steadier H-bonds, ranging from one to three, that once formed last throughout the whole simulation period. Time evolution plots of the hydrogen bonds formation and disappearance during the MD simulation along with further data including bonds identification and their occupancy can be found in Fig. S1A, B and Tables S3A–E and S4A–E.

Binding free energies, expressed as a sum of the electrostatic, Van der Waals interactions, polar and non-polar solvation free energies were calculated to better understand the compound-protein interaction during MD simulations, see Fig. S2A and B. In addition, the contribution of each protein residue to the total free binding energy was calculated. Montelukast has shown the lowest binding free energy towards 6WQF with a value of ΔGbind = -30.99 ± 5.75 kJ/mol confirming its lowest docking score among the studied compounds. The binding free energy of Ciclesonide, ΔGbind = -29.31 ± 5.35 kJ/mol, differs only slightly from the lowest value of Montelukast, with other two molecules having significantly higher (less negative) energies. However, the binding free energy of GRL-024–20, ΔGbind = -25.24 ± 5.78 kJ/mol, is slightly lower than that of Vazegepant/zevagepant, ΔGbind = -23.95 ± 5.46 kJ/mol. This is most likely caused by differences between the treatment of compound's atoms by AutoDock's and CHARMM36 forcefields. Binding free energies of complexes with PLpro (6WZU) exhibit higher fluctuations than for 3CLpro (6WQF) suggesting that these systems did not achieve equilibrium within the simulation duration (possibly due to a smaller size of the cavity), hence their averaged ensemble values are not discussed here.

Comparison with other docking studies

As previously mentioned, acquired docking results were compared to several published papers [11,[13], [14], [15], [16], [17], [18], [19], [20], [21], [22]] in order to better ascertain the compound inhibitory potential against both SARS-CoV-2 proteases. Compounds that were found in these works as well as within the already downloaded set of compounds from the PubChem database include among others well-known antivirals such as Ribavirin, Oseltamivir, Darunavir, Remdesivir or the antimalarial drug Chloroquine etc., see Table 2 . However, utilization of several different docking software with different scoring functions renders a direct comparison of binding energies ineffective.

Table 2

Docking scores of compounds found within the PubChem set and selected previously published works [11,13,[15], [16], [17], [18], [19], [20], [21]].

Identifier (CID)	Name	6LU7 ΔGbind [kcal/mol]	6WQF ΔGbind [kcal/mol]	6WZU ΔGbind [kcal/mol]	7CMD ΔGbind [kcal/mol]	Original score against 6LU7	Original score against 5R** (3CLpro) i	Original score against 3E9S (PLpro)
896	Melatonin	-7.97	-7.03	-6.02	-6.82	-10.80 a
2719	Chloroquine	-9.86	-9.06	-7.05	-5.60	-11.37 b-7.50 c-5.10 f-10.80 f		-13.88 b
6253	Cytabarine	-6.74	-6.19	-4.96	-5.79	-5.43 b		-11.70 b
37542	Ribavirin	-7.67	-6.82	-5.73	-5.23	-38.59 b		-8.30 b
44093	Captopril	-6.93	-6.10	-6.10	-5.60	-4.22 c
65028	Oseltamivir	-10.69	-9.43	-7.92	-7.42	-4.70 f	-6.891 i
73115	Clevudine	-7.02	-6.96	-5.09	-5.89	-14.52 b		15.20 b
119209	Thymidine-Methyl-T	-7.26	-6.93	-6.22	-5.52	-18.10 b		-7.57 b
131411	Arbidol	-10.40	-10.03	-6.47	-6.30	-10.12 b		-15.43 b
213039	Darunavir	-12.35	-12.92	-8.85	-6.02	-15.70 a-1.04 b-7.49 c		-8.23 b
439153	Dihydrocozymase	-12.60	-10.79	-8.69	-3.40	-11.016 d-7.0415 e
464205	Tenofovir	-9.05	-7.07	-6.06	-4.69	-25.06 b	-6.692 i	-8.49 b
492405	Favipiravir	-5.18	-4.81	-4.88	-4.41	-19.91 b-5.40 f-6.90 g		-11.89 b
3002977	Maraviroc	-12.27	-12.03	-8.62	-8.05	0.06 b		-12.61 b
5280804	Isoquercetin	-10.87	-11.43	-7.78	-5.18	-156.08 h
9875401	Rivaroxaban	-10.83	-9.83	-8.04	-9.16	-14.70 a
10182969	Apixaban	-10.50	-9.40	-7.22	-6.95	-16.40 a
10445549	Galidesivir	-7.46	-7.72	-5.87	-4.47		-6.861 i
44205240	Baricitinib	-10.10	-11.13	-7.07	-7.02		-7.253 i
121304016	Remdesivir	-12.75	-11.38	-7.71	-4.41	3.68 b-5.80 c-7.215 d-6.50 f-7.60 g	-7.529 i	-7.75 b
135398641	Inosine	-7.97	-7.03	-6.02	-4.69	-10.80 a

(a) Docking scores from Fischer et al. [13].

(b) Docking scores from Wu et al. [11].

(c) Docking scores from Hosseini and Amanlou [16].

(d) Docking scores from Hall and Ji [15].

(e) Docking scores from Owis et al. [19].

(f) Docking scores from Narkhede et al. [18].

(g) Docking scores from Silva Arouche et al.[20].

(h) Docking scores from Adem et al. [17].

(i) Docking scores from Shah et al. [21], expressed as an average of up to five docking scores against 3CLpro structures (5R7Y, 5R7Z, 5R80, 5R81, 5R82).

Therefore, 45 various best scoring compounds have been selected from the five different publications mentioned [11,13,16,17,21] and redocked using our docking protocols against all three protein targets (original selection of 52 best scoring compounds contained six that have already been evaluated in our set and one compound has been found present in two publications), see Table 3 . Several of these compounds, especially from papers [11,13,21], have indeed the calculated binding energy below -13 kcal/mol in the case of 3CLpro and below -9 kcal/mol in the case of PLpro while having a reasonable low molar mass making them well suited for a medical use, see Table 4. Candidates with the highest potential to bind with 3CLpro protease include four antiretroviral protease inhibitors, Nelfinavir [80], Lopinavir [81] and Indinavir [82]. Human CYP3A protein inhibitor Cobicistat [83], as well as the antivirals Telaprevir [84] and Simeprevir [85], used in the treatment of hepatitis C virus, have also achieved excellent docking scores. Furthermore, several flavonoids, such as Rutin, Apiin and Diosmin, and an anti-neoplastic drug Docetaxel, whose application include treatment of breast, head and neck, stomach, prostate and non-small-cell lung cancers [86] all achieved scores below -14 kcal/mol against either 6WQF or 6LU7 3CLpro structure. The highest affinity towards PLpro 6WZU structure was observed in the case of Suramine, with docking score as low as -18.92 kcal/mol, with gap of almost 8 kcal/mol to next best scoring compound. However, pose of the compound within the active site of the protease suggests cyclisation of the molecule itself, see Fig. S3, which may be the cause of the unusually low docking score and this result should be taken with a grain of salt. High affinity towards PLpro 6WZU structure was observed in the cases of abovementioned Stavudine Triphosphate, Apiin, Docetaxel and Ferristatin II. Compounds from the other studies with lowest docking scores against PLpro 7CMD structure include Diacetylcurcumin, Alectinib, Paritaprevir and two compounds without trivial names, CIDs: 101,502,236 and 101,526,067. Suramine's tendency to form intramolecular interactions was observed in smaller degree for 7CMD protein target as well (-13.74 kcal/mol). Even though, the best scoring compounds for PLpro structures form two distinct sets, the differences in their scores are well within the AutoDock scoring function error and they should be looked upon as one large set. These results have been deposited into previously mentioned excel documents containing docking results of the PubChem set of compounds and their graphical representations can be found in Fig. 7 A–D.

Table 3

Docking scores of redocked best scoring compounds from selected publications [11,13,16,17,21].

Identifier	Trivial name	6WQF ΔGbind [kcal/mol]	6LU7 ΔGbind [kcal/mol]	6WZU ΔGbind [kcal/mol]	7CMD ΔGbind [kcal/mol]	Reference
146131763	- b	-9.41	-8.41	-7.39	-7.89	[11]
CP02 a	- b	-10.54	-9.06	-7.52	-8.29
CP03 a	- b	-11.31	-9.62	-8.13	-9.16
5,750,099	- b	-11.67	-10.06	-8.15	-8.80
4,666,447	- b	-12.09	-11.47	-7.31	-7.92
CP06 a	- b	-11.50	-10.89	-7.73	-8.29
CP07 a	- b	-8.96	-8.03	-6.84	-9.16
CP08 a	- b	-11.36	-10.00	-7.89	-8.44
CP09 a	- b	-9.13	-8.17	-6.78	-8.46
CP10 a	- b	-10.06	-9.18	-7.38	-8.95
CP11 a	- b	-10.85	-9.08	-7.45	-8.38
439533	(-)-Taxifolin	-9.34	-8.75	-6.35	-7.44
4463	Nevirapine	-7.91	-7.25	-6.17	-6.26	[13]
5361	Suramine	-13.36	-15.53	-18.92	-13.74
64143	Nelfinavir	-13.02	-14.43	-7.67	-9.72
64993	6-Hept	-9.64	-8.96	-7.73	-7.42
65015	Plerixafor	-10.64	-11.23	-6.85	-7.29
65355	Stavudine Triphosphate	-9.73	-8.38	-11.40	-9.25
3010818	Telaprevir	-14.08	-13.44	-9.99	-8.90
25151504	Cobicistat	-13.54	-14.66	-9.93	-8.94
49773361	Beclabuvir	-12.71	-12.85	-7.39	-8.91
35314	O-Methyl-Dcba	-7.40	-6.86	-6.88	-6.98	[16]
107771	Pyronaridine	-13.92	-13.60	-8.99	-9.89
148124	Docetaxel	-15.12	-14.00	-10.39	-8.36
5330286	Palbociclib	-12.36	-12.23	-9.29	-9.47
9854073	Cabazitaxel	-14.85	-14.87	-9.77	-9.76
24873435	Simeprevir c	-15.70	-12.67	-9.88	-8.88
49806720	Alectinib	-11.66	-10.42	-8.24	-10.18
10621	Hesperidin	-14.04	-12.66	-9.38	-8.76	[17]
5280746	Apiin	-14.15	-13.28	-10.41	-9.01
5280805	Rutin	-14.43	-13.71	-9.26	-8.67
5281613	Diosmin	-13.88	-13.48	-9.00	-8.94
5281672	Myricetin	-11.32	-10.17	-7.31	-8.44
6441419	Diacetylcurcumin	-12.92	-12.30	-9.98	-11.38
101502236	- b	-12.38	-11.39	-9.68	-10.53
101526067	- b	-12.68	-12.12	-9.67	-10.54
92727	Lopinavir	-14.26	-15.02	-8.81	-9.37	[21]
123794	CGP42112A	-11.95	-10.25	-7.28	-7.89
392622	Ritonavir	-12.76	-13.17	-9.35	-9.81
667492	Marboran/Methisazone	-7.97	-7.27	-7.77	-7.39
5362440	Indinavir	-15.66	-14.78	-9.75	-8.65
5702678	Ferristatin II	-11.74	-10.68	-10.25	-9.25
45110509	Paritaprevir	-14.36	-14.27	-9.17	-10.28
71661251	Elbasvir	-12.42	-9.26	-7.05	-8.63
137967122	Asunaprevir	-11.99	-12.5	-7.61	-9.47

(a) Compound does not have a CID identifier.

(b) Compound does not have a trivial name.

(c) Simeprevir has been found in two different publications [16,21].

Fig. 7

Docking scores of initially selected compounds (colored diamonds), redocked compounds from selected publications (empty triangles) [13], (full triangles) [11], (empty squares) [16], (full squares) [17] and (empty circles) [21] against 6WQF (A), 6LU7 (B), 6WZU (C) and 7CMD (D) with respect to their molecular weight. Compounds with molecular weight exceeding 1000 Da, one from Wu et al. [11] (CID: 5361) and one from Shah et al. [21] (CID: 123794), are excluded from these representations.

Discussion

Highest affinity towards 3CLpro from PubChem database [26,27] was observed for anti-asthmatic drug Montelukast (CID: 5281040) [69] (6WQF) and for anti-inflammatory drug Ebastine (CID: 3191) [72] (6LU7). Interestingly, while Montelukast has reached top ten against both 3CLpro structures, Ebastine achieved a score of only -12.99 kcal/mol, ranking 18th against 6WQF structure. The anti-inflammatory immunosuppressant Solumedrol (CID: 16923) achieved the lowest docking score against PLpro 6WZU structure but was dislodged from the protein during subsequent MD simulations. Other compounds with low docking scores against both 3CLpro structures include potential anti-migraine drug Vazegepant/zevagepant (CID: 53472683) [68] and the compound GRL-024–20 (CID: 154573806).

The novel prodrug TAK-599 (CID: 73425380) achieved low docking scores against PLpro indicating its’ possible use in the COVID19 therapy. Docking protocols also revealed one compound (CID: 122146) with high affinity towards PLpro and 3CLpro 6LU7 structure that medical use has yet to be examined experimentally. Relatively high affinity towards both protein targets was also observed for Vazegepant/zavegepant, indicating that these compounds could be used as potential multi-targeting drugs. It is worth noting that many drugs that were discussed as possible SARS-CoV-2 Mpro inhibitors, such as Remdesivir, chloroquine and hydroxy-chloroquine have not achieved outstanding docking scores [11,15,16,[18], [19], [20]].

It is important to note, that most of best scoring compounds do not satisfy selected drug likeness empirical rules due to their molar mass, number of hydrogen bond acceptors/donors and/or molar refractivity. However, it is also equally important to note, that these compounds are already approved drugs that were originally proposed for other treatments, thus constituting viable set for drug refurbishment.

Comparison of obtained data with previously published results serves as an illustration of their robustness. Best scoring compounds against 6WQF show the same level of affinity towards the target as the lowest docking score compounds from other publications. Our results are on par with best results from Adem et al. [17] and have comparable molecular weights. On the other hand, redocking these structures against 6LU7 shows a decrease in docking scores compared to our set of compounds but comes at a prize of molecular weights exceeding 700 Da. Admittedly, comparison of compounds with the same molecular weight indicate that chosen set of structures obtained from the PubChem database are lagging behind best scoring compounds from other publications [11,13,[15], [16], [17], [18], [19], [20], [21]] by about 0.5 kcal/mol against 6LU7. Results of redockings against 6WZU show that compounds from the PubChem database achieve comparable docking scores with considerably lower molecular weights making them interesting targets for potential lead-like structures. However, interested reader is advised to judge the compounds not only by their docking score but also by the ligand efficacy. This parameter, expressed as the docking score with respect to the number of non-hydrogen atoms of the compound, attempts to eliminate the shift of docking score towards the more negative values caused by decrease in the Van der Waals interaction energy that goes hand in hand with increase of compound's molecular weight. Upon closer examination of ligand efficacies of the redocked compounds it becomes clear that several of the best scoring compounds are biased towards lower docking scores due to their molecular weights, see Fig. S4A–D.

In addition, we provide a set of previously published redocked compounds with reported half maximal inhibitory concentration (IC50) values against 3CLpro [38,73,87,88] and PLpro [40,[89], [90], [91], [92]] to give the interested reader more detailed view on the docking score vs. dose concentration relation of in vitro tested compounds. Refer to the supplementary materials in respective .xlsx documents, sheets IC50_dockings. However, no direct comparison has been made due to many challenges regarding such relation, e.g. imperfection of docking scores (derived from desolvation and entropy terms in AD scoring function [93]), hand in hand with relatively high uncertainties of experimental IC50 values.

Conclusion

We have performed in silico docking studies of over 860 compounds from external PubChem database [26,27] against SARS-CoV-2 proteases. Conformation stability of complexes formed by compounds with the highest affinity towards inhibition targets was validated by MD simulations and by the estimation of free energy of binding using the MM-PBSA approach [59,60]. Acquired results were compared with several previously published works [11,13,[15], [16], [17], [18], [19], [20], [21]] and show that top scoring structures are comparable with hand-picked selection of known antivirals. Docking scores database, along with ADME parameters of studied compounds, can serve for further analysis of the antiviral activity and/or for molecular design of potential drug candidates. Furthermore, the submitted data could serve as a source of validation for other researches in the fight against COVID19 and call for further in vitro verifications of antiviral activity of several best scoring compounds.

Additional supporting information

Supporting_info_README.pdf – detailed information about the supporting information files, including description of individual columns in the excel files.

Supporting_info.pdf – text document containing additional figures and tables referred to in this work.

Docking_6WQF/6LU7/6WZU/7CMD.xlsx – excel documents containing the docking results computed in this work.

PubChem_SwissADME.xlsx – excel document containing physicochemical descriptors and ADME parameters of the studied compounds.

CRediT authorship contribution statement

Marek Štekláč: Methodology, Visualization, Formal analysis, Writing – review & editing. Dávid Zajaček: Investigation, Visualization, Writing – original draft. Lukáš Bučinský: Conceptualization, Resources, Software, Writing – original draft.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.