Chemical-informatics approach to COVID-19 drug discovery: Monte Carlo based QSAR, virtual screening and molecular docking study of some in-house molecules as papain-like protease (PLpro) inhibitors.

World Health Organization characterized novel coronavirus disease (COVID-19), caused by severe acute respiratory syndrome (SARS) coronavirus-2 (SARS-CoV-2) as world pandemic. This infection has been spreading alarmingly by causing huge social and economic disruption. In order to response quickly, the inhibitors already designed against different targets of previous human coronavirus infections will be a great starting point for anti-SARS-CoV-2 inhibitors. In this study, our approach integrates different ligand based drug design strategies of some in-house chemicals. The study design was composed of some major aspects: (a) classification QSAR based data mining of diverse SARS-CoV papain-like protease (PLpro) inhibitors, (b) QSAR based virtual screening (VS) to identify in-house molecules that could be effective against putative target SARS-CoV PLpro and (c) finally validation of hits through receptor-ligand interaction analysis. This approach could be used to aid in the process of COVID-19 drug discovery. It will introduce key concepts, set the stage for QSAR based screening of active molecules against putative SARS-CoV-2 PLpro enzyme. Moreover, the QSAR models reported here would be of further use to screen large database. This study will assume that the reader is approaching the field of QSAR and molecular docking based drug discovery against SARS-CoV-2 PLpro with little prior knowledge.Communicated by Ramaswamy H. Sarma.

Introduction

Ligand based and target based drug discoveries are the most applied modern drug discovery outlooks (Hajjo et al., 2012; Polishchuk, 2017). Traditionally medicinal chemists wielded “chemical intuition” for lead optimization, which sometimes rendered biased fingerprint/fragment or scaffold selection (Bandyopadhyay et al., 2019). Now-a-days the comparative learning between the statistical and intelligent approaches enriched lead optimization aspects by combining statistical significance and putative visualizations of fingerprint (or scaffold) devoid of traditional selection-biases (Adhikari et al., 2016; Jain et al., 2020). Significantly, medicinal chemists put their primary focal point on virtual compound libraries rather than industry chemical collections to trigger time- and money-efficiency (Choudhury, 2020; Van Hilten et al., 2019).

The use of such drug discovery approaches in terms of quantitative structure-activity relationship (QSAR), artificial intelligence (AI), virtual screening (VS), drug repurposing etc. demands more when the world faced unwanted and uncontrolled scenario as like current pandemic posted by novel coronavius (2019-nCoV) (Adeoye et al., 2020; Bhardwaj et al., 2020; Elasnaoui & Chawki, 2020; Elfiky, 2020a; Mittal et al., 2020; Paniri et al., 2020; Pant et al., 2020; Patil et al., 2020; Sarma et al., 2020; Wahedi et al., 2020). In this easily accessible world of 21st century coronavirus disease-2019 (COVID-19) has been spreading alarmingly by causing huge social and economic disruption (Aanouz et al., 2020; Arya & Dwivedi, 2020; Basit et al., 2020; Boopathi et al., 2020; Ghosh et al., 2020; Hendaus, 2020; Hendaus & Jomha, 2020; Pillaiyar et al., 2020). COVID-19 respiratory disease is attributed as world pandemic by World Health Organization (https://www.who.int/dg/speeches/detail/who-director-general-s-opening-remarks-at-the-media-briefing-on-covid-19—11-march-2020). The 2019-nCoV infection has spread over to 216 countries and territories since its outbreak in the last month of 2019 in China, so far 6 287 771 confirmed cases and 379 941 deaths have been documented as on 3rd June 2020 (https://www.who.int/emergencies/diseases/novel-coronavirus-2019).

2019-nCoV is 3rd human coronavirus (HCoV) as identified in the 21st century (Ahmad et al., 2020; Anwar et al., 2020; Borkotoky & Banerjee, 2020; Chandra et al., 2020; Pillaiyar et al., 2020). Previously, two corona virus diseases severe acute respiratory syndrome (SARS) coronavirus (SARS-CoV) in 2002 and Middle East respiratory syndrome (MERS) coronavirus (MERS-CoV) in 2012 infected at least 8,422 people (fatality rate of ∼10%) and 1,700 people (fatality rate of ∼36%), respectively (http://www.who.int/csr/sars/archive/2003_05_07a/en).

2019-nCoV (also known as SARS-CoV-2; we will use SARS-CoV-2 in rest of the paper) is an envelope virus having a single-stranded positive sense RNA genome (Elfiky & Azzam, 2020; Ghosh et al., 2020; Nejadi Babadaei et al., 2020a, Babadaei et al., 2020b). The spike protein (S) of SARS-CoV-2 regulates viral entry into the host cells (Amin & Abbas, 2020; AP and VS 2020; de Oliveira et al., 2020; Elfiky, 2020b; Enayatkhani et al., 2020; Gupta et al., 2020; Hasan et al., 2020; Sinha et al., 2020; Sk et al., 2020). Two polyproteins i.e. pp1a and pp1ab are promptly translated upon entry into the host cells. Then these are disbanded by two viral proteases, one is 3 C-like protease (3CLpro) and another is papain-like protease (PLpro) (Figure 1) enzymes (Báez-Santos et al., 2015; Elmezayen et al., 2020; Ghosh et al., 2020; Joshi et al., 2020; Khan, Ali, et al., 2020; Khan, Zia, et al., 2020; Lin et al., 2018; Muralidharan et al., 2020). Both proteases are essential for SARS-CoV-2 viral replication and thus, can be is considered as druggable targets (Ghosh et al., 2020).

Figure 1.

Schematic representations of current work designwhich are composed of two major aspects- (A) Ligand based and (B) Structure based approaches.

The molecular docking and target based virtual screening studies have moved at a much faster pace (Al-Khafaji et al., 2020; Das et al., 2020; Enmozhi et al., 2020; Gyebi et al., 2020; Islam et al., 2020; Khan, Jha, et al., 2020; Kumar et al., 2020; Lobo-Galo et al., 2020; Mahanta et al., 2020) after deliberation of the first ligand bound SARS-CoV-2 3CLpro crystal structure in February, 2020. However, SARS-CoV-2 PLpro ligand based as well as structure based screening approaching were limited due to little proteomic knowledge.

In response against the social and economic disruptionposted by SARS-CoV-2 outbreak, screening of SARS-CoV PLpro inhibitors is fastest options which offer more strategic and economic benefits. The use of virtual compound libraries already gained in publicity and has achieved some successes. Here, the virtual screening has fostered the application of drug design to the SARS-CoV-2 targets. This current communication, a component of our rational drug design and discovery headway (Adhikari et al., 2017; Amin et al., 2018; Banerjee et al., 2020; Dutta et al., 2019; Halder et al., 2013), we propounded mathematical modelling workflow based on Monte Carlo optimization and other approaches which further leads to the screening of possible SARS-CoV-2 PLpro inhibitors (Figure 1).

The study design was composed of two major aspects-(A) Ligand based approaches: (i) classification QSAR based data mining of diverse SARS-CoV papin-like protease (PLpro) inhibitors, (ii) QSAR based virtual screening (VS) to identify in-house molecules that could be effective against putative target SARS-CoV PLpro; and (B) Structure based approaches: finally validation of hits through receptor-ligand interaction analysis (Figure 1).

Therefore, this study may introduce key concepts, set the stage for molecule identification and QSAR based screening of in-house molecules active against putative SARS-CoV-2 PLpro enzyme.

Methods and materials

Dataset

A set of diverse SARS-CoV PLpro inhibitors were collected with inhibitory activities (Báez-Santos et al., 2014; Cheng et al., 2015; Chou et al., 2008; Frieman et al., 2011; Ghosh et al., 2009, 2010; Park et al., 2012; Ratia et al., 2006, 2008). Compounds with no inhibitory activity and without definite activity were not taken for this study. In addition, duplicate molecules were eliminated. Finally, ninety one molecules were considered for the further molecular modelling study (Table S1). The average SARS-CoV PLpro pIC value was considered as the ‘activity threshold’ for the current study. Compounds having the SARS-CoV PLpro pIC value less than the ‘activity threshold’ were classified as lower PLpro inhibitors or Inactives and those with PLpro pIC value higher than the ‘activity threshold’ were yielded as promising PLpro inhibitors or Actives. Thus, out of 91 molecules, 40 compounds were distinguished as Actives and 51 molecules were considered as Inactives.

Classification based QSAR

The classification modelling assists to discriminate the Active and Inactive molecules in terms of their investigated biological significance. Here, we performed Monte Carlo based Coral QSAR study. Performing this study not only offers a graphical visualization of critical fingerprint or fragments attributed to enhance/decrease the SARS-CoV PLpro inhibitory activity but also it allows the chance of screening external set compounds.

Monte Carlo optimization based QSAR

Descriptors calculation

SMILES-based descriptors

SMILES-based descriptors were calculated by the following equation:

In this equation, T and N represent threshold value and number of epoch, respectively. The correlation weights were represented by CW. The different coefficients like a, b, c, d, α, β and γ were used for descriptor modification. NOSP, HALO, BOND and ATOMPAIR represent global SMILES attributes and the local smile attributes were denoted by Sk, SSk and SSSk (Toropov et al., 2013; Toropova et al., 2015).

Graph-based descriptors

GAO (graph of atomic orbital), HSG (hydrogen-suppressed graph) and HFG (hydrogen-filled graph) represents different graph based descriptors and was calculated by following equation:

Where, 0ECk, 1ECk and 3ECk represent different Morgan’s connectivity indices. Ak denotes different chemical atoms like: C, N, O etc. α, β and γ were the coefficients with 0 and 1 value. The coefficients having value 0 and 1 were denoted by as α, β and γ (Toropov et al., 2013; Worachartcheewan et al., 2014).

Hybrid descriptors

The amalgamation of SMILES and graph-based descriptors forms hybrid descriptors which are represented as:

Model development and validation

In our study, by using balance of correlation method twenty-one classification models were developed from three different splits. The dataset containing of 91PLpro inhibitors were distributed into training (41 compounds), calibration (35 compounds) and test (15 compounds) sets which were used for the study. Further, optimization of T (threshold) and N (epoch) were also performed separately for each model (Toropova et al., 2015). The sensitivity, specificity, accuracy along with the MCC values was recorded as a measure of internal and external validation. Finally, the important structural attributes that were solely answerable for promoting or hindering of PLpro activity were identified.

Target based molecular modelling

Homology modelling was performed which provided 3 D models for SARS-CoV-2 protein structure as the ligand-bound crystal structures are not available till date. The homology model for SARS-CoV-2 was built using Swiss Model web server (https://swissmodel.expasy.org/) and subsequently, validated by Verify3D (https://servicesn.mbi.ucla.edu/Verify3D/), ProSA (https://prosa.services.came.sbg.ac.at/prosa.php) and PROCHECK (https://spdbv.vital-it.ch/). In addition, the model was optimized using Swiss PDB viewer software using GROMOS96 Force-Field followed by determination of RMSD value by the aid of PyMOL software (https://pymol.org/2/). Lastly, the energy minimized model was implemented for the molecular docking analysis. The docking study was conducted by the aid of AutoDock Vina (Trott & Olson, 2009). Notably, grid box was selected by covering the geometric pattern occupied by the prototype in-bound ligand in the crystal structure of PDB: 4OW0. The docked poses of ligands were visualized by PyMOL software (https://pymol.org/2/) and the 2 D-interaction plots were generated by Discovery Studio 3.5 Visualizer (Accelrys Software Inc., San Diego, California, USA).

Result and discussions

It is already reported that some small molecules exhibited potent SARS-CoV PLpro inhibition (Ghosh et al., 2009, 2010, 2020).SARS-CoV PLpro shares 82.80% sequence similarity with the homologus SARS-CoV-2 strain (Figure 2).

Figure 2.

Comparison of SARS-CoV PLpro and SARS-CoV2 PLpro: A modeled structure of SARS‐CoV‐2 PLpro (orange), the crystal structure of SARS‐CoV (grey, PDB: 4OW0) where PLpro inhibitor binds at the catalytic site [docked ligand (brick red) and in-bound ligand (light yellow)].

Significantly, PLpro active site amino acids (P248, P249, Y269, D165, E168, L163, G164, Q270, Y274, Y265, T302) of both stains are highly conserved (Figure 2). Thus, it may pretend that SARS-CoV PLpro inhibitors would to be potential inhibitors of SARS-CoV-2 PLpro enzyme. As the development of new small molecules against the proteasesof the COVID-19 is challenging as well as time and money consuming, it is better to screen compounds based on the previous ones.

The Swiss model constructed an excellent homology model of SARS-CoV-2 based on the sequence identity between the PLpro SARS-CoV (PDB: 4OW0) and PLpro SARS-CoV-2. The quality of the model was validated by the aid of Verify3D (96.79% of the residues had average 3 D -1 D score s ≤ 0.2), ProSA server (Z-Score = −8.79), the Ramachandran plot (91.4% and 8.2% of the residues in the most favorable and the additional allowed region, respectively; 0.4% residues in generously allowed regions while no residues were found to be in disallowed region). The RMSD score, as determined by using PyMOL, was recorded 0.090 which confirmed the model acceptability.

The amino acid sequences of SARS-CoV PLpro (PDB: 4OW0) and COVID-19 PLpro (homology modelled) are depicted. Notably, the active site amino acids including P248, P249, Y269, D165, E168, L163, G164, Q270, Y274, Y265, T302 etc. are highly conserved.

Classification QSAR study

The whole set of molecules able to bind to the SARS-CoV PLpro enzyme were taken after extensive literature studies, retrieving only those ligands having an absolute IC50 values. This set consisted of 91 PLpro inhibitors. Depending on the ‘activity threshold’, out of 91 compounds, 40 compounds were identified as Active and 51 molecules were considered as Inactives.

Monte Carlo optimization based coral QSAR

In Monte Carlo optimization (Toropov et al., 2018; Toropova et al., 2015), a total of twenty-one different models from three different splits were generated using SMILES and graph-based descriptors with a combination of different connectivity indices which were calculated for generation of different models (Table S2). Each models were developed after the search for desirable T (threshold) and N (epoch) values as per the test set statistics as suggested by Toropova et al. (2015). The statistical parameters of three best models from three different splits are shown in Table 1.

Table 1.

The statistical performance of three best models from three different splits.

Parameter	Set	T	N	TP	TN	FP	FN	Total	Sensitivity	Specificity	Accuracy	MCC
M2 SMILES, GAO (0ECk)	Sub-Training	3	7	15	24	2	0	41	1.0000	0.9231	0.9512	0.9025
	Calibration			15	18	0	2	35	0.8824	0.1000	0.9429	0.8911
	Test			4	4	3	4	15	0.5000	0.5741	0.5333	0.5014
M13 SMILES, HSG (0ECk)	Sub-Training	1	7	15	24	0	2	41	0.8824	1.0000	0.9512	0.9025
	Calibration			17	17	1	0	35	1.0000	0.9444	0.9714	0.9444
	Test			5	7	2	1	15	0.8333	0.7778	0.8000	0.6001
M15 SMILES	Sub-Training	2	6	10	30	0	1	41	0.9091	1.0000	0.9756	0.9380
	Calibration			20	15	0	0	35	1.0000	1.0000	1.0000	1.0000
	Test			8	4	2	1	15	0.8889	0.6667	0.8000	0.5774

Where, T = Threshold, N = Epoch;The selected model is shown in bold face,.

Among all the 21 models developed in 3 different splits, model M13 (SMILES and HSG with 0ECk) from split-2 was found to be the best one.

The end point values calculate for M13 is shown below:.

Endpoint = 0.1282 (± 0.0048) + 0.04293 (± 0.00023) * DCW (1,7).

As can be seen from Table 1, the model M13 showed a satisfactory predictive ability. The values of the sensitivity, specificity, accuracy along with the MCC obtained for both the sub-training and calibration sets were highly encouraging. Indeed, the values attained for the test set (i.e. sensitivity, specificity and accuracy of 0.8333, 0.7778 and 0.8000, respectively) suggested the acceptable external predictive power of the classification based model. However, the MCC value of the test set was comparatively poor than the MCC values of the sub-training and calibration sets. Different structural attributes of the best model M13 (SMILES and HSG with 0ECk) from split-2 is depicted in Table S3.

QSAR based virtual screening

Since the QSAR models proved reasonably accurate, we used the best model to screen our in-house molecules (IH-001 – IH-067) from our previous publications (Adhikari et al., 2016; Halder et al., 2015; Mukherjee et al., 2017). General structures of our in-house compounds (IH-001 – IH-067) are depicted in Figure 3.

Figure 3.

General structures of our in-house compounds (IH-001 – IH-067).

These compounds were already reported as metalloproteinase inhibitors (Adhikari et al., 2016; Halder et al., 2015; Mukherjee et al., 2017). First, we predicted the compounds and then screening as per their potentiality in the Monte Carlo based classification QSAR model. Second, we defined an applicability criterion to choose the best hits.

In this regards, we first screened the sixty seven in-house compounds. From the in-house database, a collection of 56 compounds were predicted as Active from the Monte Carlo based QSAR model (Table S4). After screening 56 similar compounds in SwissADME (Daina et al., 2017) − 13 compounds including IH-009, IH-015, IH-017, IH-020, IH-023, IH-027, IH-037, IH-038, IH-040, IH-043, IH-046, IH-047 were passed the ADME criteria (Figure 4).

Figure 4.

Radar plot of the in-house compounds after calculating ADME data by SwissADME server (http://www.swissadme.ch/) suggesting the drug-likeness [the pink area represents the optimal range of each properties.LIPO = Lipophilicity (between −0.7 and +5.0), SIZE = Molecular weight (between 150 and 500 g/mol), POLAR = Polarity (between 20 and 130Å2), INSOLU = Solubility (not higher than 6), INSATU = Saturation (fraction of carbons in the sp3 hybridization not less than 0.25), FLEX = Flexibility (no more than 9 rotatable bonds)].

Lastly, these potential chemotypes were considered for molecular docking study against the putative target SARS-CoV-2 PLpro.

Binding interactions of the in-house VS hits to SARS-CoV-2 PLpro

Molecular docking study offers a vital tool to predict the possible structural conformations between ligand and active sites of a receptor/enzyme. Here, molecular docking approach was employed on the in-house VS hits by using the AutoDock Vina to understand the docking/binding mode between in-house VS hits and SARS-CoV-2 PLpro. Meanwhile, the docking modes of two prototype in-house VS hits in the catalytic site amino acid residues of COVID-19 PLpro are illustrated in Figure 5.

Figure 5.

The docking modes of two prototype in-house VS hits (A) IH-009 and (B) IH-027 in the catalytic site amino acid residues of SARS-CoV-2 PLpro (proetin, grey cartoon; active site amino acids, orange stick; in-bound lignad, red stick; in-house molecules, yellow stick; water molecule, red ball; hydrogen bond interactions, light green dashed lines).

Furthermore, the re-dock binding pose of the in-bound PLpro inhibitor (red stick) on the SARS-CoV-2 active site are shown in Figure 5. Notably, the positions of the in-house VS hits in the COVID-19 PLpro catalytic site are basically the same from Figure 5, suggesting that in-house VS hits capture the right position in the PLpro cavity. This phenomenon validated the accuracy of docking study.

The docking poses of the all thirteen in-house VS hits are depicted in Figure S1. The entire in-house VS hits snugly occupied the binding site of SARS-CoV-2 Plpro. In addition, the Figure S1 pinnacle the superimposition of docking poses of thirteen in-house VS hits in the homology modelled SARS-CoV-2 PLpro. This observation justified that these investigated derivatives exhibit potential to be PLpro inhibitor and may be a valid weapon against SARS-CoV-2.

The binding interaction of IH-009 with SARS-CoV-2 PLpro revealed three hydrogen bonds with three amino acids (Leu162, Tyr264 and Tyr268), one π-π T-shaped interaction between phenyl ring and Tyr268, along with additional π-alkyl interaction between phenyl ring and Pro248. Apart from that the bromine molecule of IH-009 interacted with two proline amino acid residues (Pro 247 and 248). The carboxylic acid feature formed water mediated hydrogen bond interactions. These interactions were more or less consistence with other molecules also (Table S5). Notably, the napthyl ring formed two π-π T-shaped and two π-alkyl interactions. An additional π-sigma interaction was noticed where Tyr264 was involved as given in Table S5. However, the exact mechanism of the binding is still sketchy as it required further molecular dynamic simulation study. Moreover, the in vivo effects of these in-house VS hits would be needed to confirm the mechanism.

Conclusion

Here, we have constructed a classification based QSAR model that could be used as a tool for predicting new molecules and/or virtual screening. The model developed by Monte Carlo optimization based QSAR were followed by virtual screening of some in-house chemicals. Then ADME data driven screening was performed by SwissADME and identified compounds with good drug-likeness. Finally, molecular docking study of QSAR derived virtual hits was performed to increase the confidence in the final hypotheses. The molecular docking study performed against putative target SARS-CoV-2 PLpro suggesting potentiality of these investigated in-house molecules. Thus, it can be concluded that the in-house molecules have potential to use as a seed for drug design and optimization against SARS-CoV-2 PLpro. After extensive in vitro and in vivo studies, these in-house VS hits may be emerged as therapeutic options for COVID-19. This study may also motivate medicinal chemists to design similar type of compounds in hopes to trigger biological potency as well as efficacy without accruing toxicities.