key: cord-0952323-kgtx9pzf
authors: Patel, Dhaval; Athar, Mohd; Jha, P. C
title: Computational investigation of binding of chloroquinone and hydroxychloroquinone against PLPro of SARS-CoV-2
date: 2020-11-17
journal: Journal of biomolecular structure & dynamics
DOI: 10.1080/07391102.2020.1844804
sha: 51ba7ecdd6fdc24f369bff1e14846e08a2b337e0
doc_id: 952323
cord_uid: kgtx9pzf

Novel coronavirus SARS-CoV-2 has infected 18 million people with 700,000+ mortalities worldwide and this deadly numeric figure is rapidly rising. With very few success stories, the therapeutic targeting of this epidemic has been mainly attributed to main protease (Mpro), whilst Papain-like proteases (PLpro) also plays a vital role in the processing of replicase polyprotein. Multifunctional roles of PLpro such as viral polypeptide cleavage, de-ISGlyation and immune suppression have made it a promising drug target for therapeutic interventions. Whilst there have been a number of studies and others are on-going on repurposing and new-small molecule screening, albeit previously FDA approved drugs viz. Chloroquine (CQ) and Hydroxychloroquine (HCQ) have only been found effective against this pandemic. Inspired by this fact, we have carried out molecular docking and dynamics simulation studies of FDA approved CQ and HCQ against SARS-CoV-2 PLpro. The end aim is to characterise the binding mode of CQ and HCQ and identify the key amino acid residues involved in the mechanism of action. Further, molecular dynamics simulations (MDS) were carried out with the docked complex to search for the conformational space and for understanding the integrity of binding mode. We showed that the CQ and HCQ can bind with better binding affinity with PLpro as compared to reference known PLpro inhibitor. Based on the presented findings, it can be anticipated that the SARS-CoV-2 PLpro may act as molecular target of CQ and HCQ, and can be projected for further exploration to design potent inhibitors of SARS-CoV-2 PLpro in the near future.

The year 2019 ended with a fatal outbreak of a novel coronavirus (SARS-CoV-2) preliminary identified as a causative agent for a series of unusual pneumonia cases in Wuhan city, Hubei province of China (Chan et al., 2020; . As the cases rises sporadically, the WHO announced Public Health Emergency of International Concern (PHEIC) for the 2019-nCoV outbreak in January 2020 (WHO.int). Since the infection crossed geographical barriers, the WHO permanently named the 2019-nCoV pathogen as SARS-CoV-2 (coronavirus disease 2019, COVID-2019) and declared pandemic situation in March 2020 (WHO.int). As of now (11 October 2020), more than 36 million confirmed cases with >1 million death are reported globally affecting >208 countries in the world (WHO.int). The disease caused by SARS-CoV-2 presents vast pathophysiological symptoms including fever, coughing and shortness of breath in common cases whereas pneumonia, kidney failure and severe acute respiratory failure in severe cases (Chan et al., 2020; .

The causative agent for SARS-CoV-2, coronavirus belongs to subgenus Sarbecovirus, Coronaviridae family and Orthocoronavirinae subfamily. It is a b-coronavirus of 2B (b-coronavirus), however the group is non-zoonotic till date unlike earlier coronaviruses Weiss & Leibowitz, 2011) . SARS-CoV-2 is an enveloped single-stranded RNA virus (þve ssRNA) with a genome size of 29.9 kb that spreads widely amongst humans and other mammals, causing a wide range of infections from common cold symptoms to fatal diseases, such as severe respiratory syndrome Zhu et al., 2020) . Genome sequencing of SARS-CoV-2 has revealed 96% identity to the bat coronavirus and 79.6% sequence identity to SARS-CoV . Despite being genetically distinct from SARS-CoV, the proteome is quite similar with only difference of two major proteins responsible for pathogenicity: the non-structural proteins (nsps) and structural proteins (Naqvi et al., 2020) . With the aid of ribosomal frameshifting during translation, the replicase gene of SARS-CoV-2 encodes two overlapping translation products, polyproteins 1a and 1ab (pp1a and pp1ab) . Each of these polyproteins is then cleaved by main protease (Mpro) and Papain-like Protease (PLpro) of SARS-CoV-2 to release 11 (pp 1a) and 5 (pp 1ab) functional proteins, respectively, necessary for viral replication.

As viral proteases play a vital role in processing the polyproteins that are translated from the viral RNA, inhibitors designed against these proteases may prove effective in blocking the replications of coronaviruses, thus making them attractive target for drug discovery . Viral proteases have been suggested to be linked with the mechanism of infection and pathogenicity of SARS-CoV-2 (Benvenuto et al., 2020) . The papain-like protease (PLpro) recognises LXGG#(A/K)X sequence motif for cleavage of its polyprotein for releasing five functional non-structural proteins (nsp). Apart from its main role, PLpro has also been reported to suppress host innate immune responses by subverting cellular Ubiquitination machinery (Barretto et al., 2005; Ratia et al., 2006) and also involved in de-ISGlyation for facilitating viral survival (Daczkowski et al., 2017; Lee et al., 2015) . Taken together, multifunctional role of PLpro have positioned this as target of choice for combating the viral infection (B aez-Santos et al., 2014 (B aez-Santos et al., , 2015 . Back in 2008, the first class of non-covalent drug-like naphthalene PLpro inhibitors and their analogues were discovered that exhibited nanomolar inhibition against SARS-CoV PLpro (Ghosh et al., 2009; Ratia et al., 2008) .

A number of inhibitors have been recently proposed for the treatment of COVID-19 under the 'public health emergency' situation. Amongst them, most of the studies have been focussed on drug repurposing of existing known inhibitors for previous CoVs or other viruses such as HIV, influenza and Ebola. This includes Chloroquine (CQ) and its derivative Hydroxychloroquine (HCQ) (antimalarial drugs) (Colson et al., 2020; Gao et al., 2020; Yao et al., 2020) ; Lopinavir and Ritonavir (HIV Protease anti-viral inhibitors) (Chu et al., 2004; Yao et al., 2020) as inhibitors for Mpro, Remdesivir (GS-5734) (nucleoside analogue) (Agostini et al., 2018; and Favipiravir as inhibitors of RNAdependent RNA polymerase , etc. Recent clinical studies and computational data have reported that CQ has a potential to treat COVID-19 (Adeoye et al. 2020; Boopathi et al. 2020; Gao et al., 2020) . Previously, these two drugs (CQ, HCQ) were approved by FDA for the treatment of malaria, rheumatoid arthritis, lupus and sun allergies (Boopathi et al., 2020) . The WHO have also now recommended antimalarial drug chloroquine, to treat coronavirus (WHO.int; solidarity trials). The putative anti-viral effects of CQ have been hypothesised to be related with the elevation of endosomal and lysosomal pH in addition to its interference with terminal glycosylation of SARS-CoV's S protein receptor angiotensin-converting enzyme 2 and in inhibiting sialic acid biosynthesis via quinone reductase 2 (Gay et al. 2012; Vincent et al., 2005) . The derivative of CQ, HCQ is also hypothesised to act on other viruses in a similar manner, suggesting towards the identical mechanism of action. Recent molecular docking and simulations studies have revealed the binding mechanism and stabilisation of HCQ with SARS-CoV-2 main protease (Mpro) through non-covalent interactions (Mukherjee et al., 2020) . Another similar molecular docking study of HCQ against main protease (Mpro) have suggested that compound binding is facilitated by the conformational and rotational changes (Baildya et al., 2020) . Recently, Nimgampalle & co-workers have also performed molecular docking studies of CQ, HCQ and its derivatives against multiple SAR-CoV-2 proteins (Spike glycoprotein, RNA-dependent RNA polymerase, Chimeric Receptor-binding domain, Main protease, Non-structural Protein 3, Non-structural Protein 9, ADP-ribose-1 monophosphatase) (Nimgampalle et al., 2020) . Amongst the various derivatives of CQ, CQN2H and CQN1B showed potential inhibition against all the drug targets. Studies have also suggested HCQ has more affinity compared to CQ with the NTD-N-protein of SARS-CoV-2 (Amin & Abbas, 2020) . In spite of HQ and HCQ proposition for the treatment of COVID-19, the exact molecular mechanism and target in which CQ and HCQ interact with SARS-CoV-2 is yet to be answered. Moreover, no single study has been undertaken to evaluate the binding mechanism and interactions of CQ and HCQ with PLpro of SARS-CoV-2. Owing to the need of the hour in the on-going pandemic, we have carried out molecular docking and dynamics simulation studies of FDA approved CQ and HCQ against SARS-CoV-2 PLpro. We aimed to characterise the binding mode, binding affinities and key amino acid residues playing a crucial role in their mechanism of action. Through this study, we aspired to compare CQ and HQ with known inhibitors of same class and project them for further exploration to design potent inhibitors of SARS-CoV-2 PLpro in the near future.

2.1. Sequence retrieval, homology modelling, refinement and structural validation of SARS-CoV-2 PLpro

The amino acid sequence of SARS-CoV-2 PLpro was retrieved from NCBI with GenBank accession number QJA17047.1 (305 residues) in FASTA format for building homology model using the I-TASSER server (https://zhanglab.ccmb.med.umich. edu/I-TASSER/) (Yang et al., 2015) . I-TASSER first identifies structural templates from the PDB by multiple threading approach followed by full-length atomic models construction in iterative template-based fragment assembly simulations. The model was evaluated by the Ramachandran method, as predicted with Rampage server (http://mordred.bioc.cam.ac. uk/$rapper/rampage.php) (Lovell et al., 2003) . Further, the structure was refined to improve the overall model quality by ModRefiner server (https://zhanglab.ccmb.med.umich.edu/ ModRefiner/) (Xu & Zhang, 2011 ).

The 2-D structures of the compounds were sketched using MarvinSketch program and exported in .SDF format. These ligand molecules were prepared in LigPrep module of Schrodinger suite. The LigPrep program checks for alternative tautomers, chirality, ionisation states and low energy ring conformations. The prepared ligands were considered for multiconformation generation by Macromodel conformation search module i.e. Mixed-torsional/Low-mode sampling (MTLM) at the default settings. In particular, OPLS_2005 force-field and TNCG methods was used for energy minimisation in 500 steps (Jorgensen et al., 1996; Shivakumar et al., 2010) . The Glide molecular docking of Schrodinger suite was utilised to predict the binding mode of ligands (Halgren et al., 2004) . The protein structure was processed using protein preparation wizard and the docking grid was generated using TTT ligand of PLPro (PDB: 3e9s) (Madhavi Sastry et al., 2013) . The Glide-SP (standard precision) docking protocol starts with the systematic conformational expansion of the ligand, followed by its placement in receptor site. Default parameters were used for glide docking run.

The docked complex of SARS-CoV-2 PLpro with four ligands (TTT, PB5, CQ, HCQ) were used as starting conformation in Molecular Dynamics Simulations (MDS). For assessing the binding stability and molecular interaction, apo and holo structures (docked complexes) which include CQ, HCQ and reference inhibitor (TTT and PB5) were subjected for MD simulation studies. A total of five MDS systems were prepared and subjected for MDS with 50 ns production run (SARS-CoV-2 PLpro; SARS-CoV-2 PLpro_PB5; SARS-CoV-2 PLpro_TTT; SARS-CoV-2 PLpro_CQ and SARS-CoV-2 PLpro_HCQ). MD simulation was performed with GROMACS ver. 2016 .4 (Van Der Spoel, 2005 with Amber99SB force-field (Lindorff-Larsen et al., 2010) and the protocol was essentially the same as mentioned in our previous reports (Patel et al., 2017 . Briefly, ligand parameter and topology was prepared by ACPYPE (Sousa da Silva & Vranken, 2012), wherever required. All the MDS systems were solvated with Three-site water model (TIP3P) in dodecahedron box, maintaining a distance of 1 nm from edges of protein in all directions. The MDS systems were neutralised with an equal number of counter ions (Na þ /Cl À ) followed by energy minimisation using the steepest descent algorithm to remove any steric clashes and bad contacts with maximum force <1000 kJ mol À 1 nm À 1 (50,000 steps max (Berendsen et al., 1984) was used for maintaining the system at constant volume (100 ps) and at a constant temperature (300 K). Subsequently, NPT equilibration was performed at a constant pressure (1 bar) for 100 ps maintained by Parrinello-Rahman barostat (Parrinello & Rahman, 1980) . The Particle Mesh Ewald approximation was applied with 1 nm cut-off for calculating long-range electrostatic interactions, computing coulomb & the Van der Waals interactions (Darden et al., 1993) . The bond length was constrained using the LINCS algorithm (Hess et al., 1997) . Finally, 50 ns production run was executed with default parameters, and the coordinates were saved after every 2 fs time frame. The MDS trajectories were visualised using VMD (Humphrey et al., 1996) and Chimera (Pettersen et al., 2004) . For calculation of root mean square deviation (RMSD), root mean square fluctuation (RMSF), hydrogen bonds (H-bonds) etc. 'gmx' commands were used in GROMACS and the plotting tool GRACE was used for generating the plots (http:// plasma-gate.weizmann.ac.il/Grace) . For molecular visualization and other analysis Discovery studio was utilized (BIOVIA 2020).

The dominant conformation in the entire MDS trajectory was assessed by RMSD based clustering via 'gmx cluster' which explore the conformational landscape amongst the ensemble of protein structures. The GROMOS algorithm (Daura et al., 1999) was used to determine the dominant conformation with Ca RMSD cut-off value 0.15 nm. Further, Principal Component Analysis (PCA) or Essential Dynamics (ED) analysis was carried out to understand the collective and overall motion, as reported in our previous study using 'gmx covar' and 'gmx anaeig' tools (Patel et al., 2020) . PCA reduces the complexity of the data and results in the concerted motion indicating important structural changes. The set of eigenvectors and eigenvalues were computed by diagonalising the covariance matrix after removing the translational and rotational motions. The amplitude of the eigenvector is represented by eigenvalues in multidimensional space, whilst the Ca displacement along each eigenvector shows the concerted motions of the protein along each direction. Free Energy Surface (FES) (kcal/ mol) was also computed for all the MDS systems considering the conformational variability in terms of ROG and RMSD.

Binding free energy of SARS-CoV-2 PLpro docked complexes was estimated by Poisson-Boltzmann or generalised Born surface area continuum solvation method (MM/PBSA and MM/ GBSA) (Genheden & Ryde, 2015; Sun et al., 2018; Wang et al., 2016) . Resulting energies of individual residues were used to quantitatively estimate the ligand affinity for SARS-CoV-2 PLpro. The 'g_mmpbsa' tool (Kumari et al., 2014) with default parameters was used for calculating molecular mechanics potential energy (electrostatic þ Van der Waals interactions) and solvation free energy (polar þ non-polar solvation energies) calculations. The last converged 30 ns (150 frames) trajectories of the complexes were assessed by the RMSD plot for estimating binding free energy. The frames were selected at a regular interval of 200 ps covering a wide range of conformational space.

The SARS-CoV-2 PLpro shares 83% sequence identity with the SARS-CoV PLpro and also comprise considerable percentage identity with other CoVs PLpro (Figure 1 ). Full-length homology model of SARS-CoV-2 PLpro (QKS90154) was generated using the I-TASSER server ( Figure 2 ). I-TASSER utilises multi-template library for the model construction. We observed that the C-score of the model was 1.86 whereas estimated TM score was 0.98 ± 0.05. The resulting model was refined by ModRefiner server and then subsequently evaluated for Ramachandran statistics followed by the energy minimisation. The Ramachandran plot and statistics are appended in supplementary Figure S1 . Sequence analysis suggest that there are 54 variations in the sequences of SARS-CoV-2 PLpro as compared to previous SARS-CoV PLpro ( Figure 1) . However, 40 out of 54 variations were located over the finger, palm, thumb and UbL domain of PLpro (of both CoV), as can be seen by 3-D surface mapping ( Figure 2 ). In particular, we noted one important variation in SARS-CoV-2 PLpro i.e. C104S which forms catalytic triad.

The ligands were prepared by the ligprep and their ionisation states were assigned at pH 7 ± 2 using Epik. The docking grid was generated using co-crystal PLPro ligand (TTT) taken from PDB entry 3E9S. Control docking experiments by redocking was performed and we obtained similar binding pose, subsequently; same grid was employed for docking CQ and HCQ ligand. Four possible ionisation states were generated for HCQ whereas two states for CQ at pH 7 ± 2 (and for 4ovz & 3e9s). After executing Glide docking with default parameters, top scored poses were shortlisted based on fingerprint interactions with reference ligand. Representative binding modes of the compounds has been provided in the Figure 3 .

For docking experiments, we were enthused by the fact that TTT ligand has aromatic naphthalene ring which engage in hydrophobic interactions with Tyr274, Tyr265, Thr302, Pro248, Pro248 and Pro249. Likewise, CQ and HQ also comprise quinolone ring which we believe will act as anchoring point and thereby both compounds will also pack into the similar hydrophobic pocket.

As shown in Figure 3 , PB85 (PDB: 4OVZ) can form two Hbonds, one with N of Tyr269 and another with protonated Figure 1 . Multiple structure-based sequence alignment of SARS-CoV-2 PLpro taking as a reference structure and its comparison with SARS-CoV PLpro. Secondary structure assignments are marked in the alignment output along with catalytic triad residues labelled by black stars, whilst residues forming the zinc-finger motif are marked with blue stars. NH with Asp165. TTT also interacted with Asp165 through Hbonding with NH whereas Glu168 interact with amino group of terminal benzylamine ring of TTT. The Tyr269 was commonly engaged through pi-pi interactions in both the crystal structures. Upon inspecting the crystal structures, we believe that Tyr269 orientation supports the compound binding by covering the site at the linker region of two terminal rings. This way, it prevents the compound liability for the binding. In particular, TTT makes contacts with Tyr269 which renders the wall-like support to the pocket whereby N acceptor forms H-bond with Asp165, and carbonyl donor form Hbond interaction with the side chain of Tyr265. Amino group of benzamide forms H-bond with Gln270. The detailed illustration of ligand-receptor interactions is shown in Figure 3 .

For CQ, we got the expected binding mode in which hydrophobic ring recapitulates the aromatic rings of known ligand, whereas for HCQ, we achieved distinct pose (mode 1) top scored binding mode. Nonetheless, we selected one more binding mode for HCQ (mode 2) where quinolone ring orients in the hydrophobic site. One possible reason for distinct orientation of HCQ is the presence of additional group (absent in CQ), which makes additional H-bonding interaction Arg167 in mode1 (supplementary Figure S2) , thus deprioritised the CQ-like pose. However, we tested all the poses for the molecular dynamics simulations. The docking scores and details of major interaction are given in Table 1 .

The MD simulations were carried out to investigate the stability, dynamics and conformational changes of docked PLpro-inhibitor complexes. The SARS-CoV-2 PLpro docked ligands (CQ and HCQ) and two known reference PLpro inhibitors (PB5 and TTT) were subjected for MD simulations. Later, we analysed RMSD, RMSF, H-bonds as well as PCA component of the CQ and HCQ complexes and compared with apo PLpro as well as reference ligands. Finally, the binding free energy of all complexes was computed for the last stable 30 ns trajectory. 

The stability of the PLpro-inhibitor complexes was evaluated to gain molecular insights into the binding interactions of CQ and HCQ via RMSD value of Ca backbone and ligand for the entire 50 ns simulations. The RMSD measures the structural variation between the Ca backbones from its initial conformation to its final position during the entire simulation trajectory. Thus, lesser the RMSD values indicate higher stability of the simulation system. The RMSD values for all five MD systems are shown in plot (RMSD (nm) vs. Time (ns)) in Figure 4A for SARS-CoV-2 PLpro; and its complexes with TTT, CQ and HCQ. The mean RMSD values for protein in Apo PLpro and its complexes with PB5, TTT, CQ and HCQ were 0.164 ± 0.025, 0.178 ± 0.03, 0.147 ± 0.021, 0.191 ± 0.0322 and 0.161 ± 0.025, respectively. The RMSD deviation of CQ and HCQ complexes are found to be similar in magnitude with the reference inhibitor complexes and apo protein. This connotes that the docked complexes are more stable during the simulations than the apo. As displayed in Figure 4B , the mean ligand RMSD values for PB5, TTT, CQ and HCQ were 0.123 ± 0.032, 0.145 ± 0.013, 0.099 ± 0.0299 and 0.283 ± 0.04, respectively. From the above observation, it can be concluded that CQ and HCQ behave well within the active site of the PLpro protein.

Further, the radius of gyration (ROG) was calculated for all five MDS for the entire 50 ns trajectory. ROG is a measure of protein shape at each time-point and comparable to the experimentally obtainable hydrodynamic radius. The average ROG values of protein Ca backbones for SARS-CoV-2 PLpro; and its complexes with PB5, TTT, CQ and HCQ were 2.27 ± 0.016, 2.26 ± 0.016, 2.26 ± 0.016, 2.26 ± 0.018 and 2.26 ± 0.018 nm, respectively (supplementary Figure S5) . The above ROG values of all the MDS systems were similar indicating that the overall shape of the apo and inhibitor complex PLpro are consistent irrespective of the inhibitor.

Fluctuations in residues evolved in 50 ns trajectory were computed from the Root Mean Squared Fluctuations (RMSF) plot ( Figure 5A ). The RMSF values of each PLpro-inhibitor complexes were overlaid on the apo PLpro for comparison of flexible residues in the absence and presence of inhibitor ( Figure 5A ). The mean RMSF value for Apo_PLpro; PB5; TTT; CQ and HCQ were 0.123 ± 0.066, 0.125 ± 0.066, 0.123 ± 0.057, 0.127 ± 0.071 and 0.123 ± 0.054 respectively. Overall, the RMSF values of PLpro in complex with CQ and HCQ were very similar in the pattern when compared to apo PLpro and PLpro complex with reference ligand. Though a very few fluctuations were observed at some residues. In particular, the residues in range 190-200 and 225-230 showed higher fluctuations in case of apo PLpro and CQ-PLpro complex as compared to HCQ-PLpro complex and reference ligand complexes (PB5 and TTT). Whereas, the residue at position 270 showed high fluctuations for apo PLpro as compared to PLpro complexes with inhibitors, likewise, a similar fluctuations pattern was also found for TTT and HCQ. The residues involved in catalytic triad formation (S104, H264 and D279) and residues forming zinc-finger motif (C182, C185, C217 and C219) had identical RMSF values indicating the conserved stability in absence and presence of inhibitors. The ligand RMSF values were calculated to assess the overall fluctuations of the ligand atoms ( Figure 5B ). These fluctuations were in a similar range as compared with reference ligand.

The RMSF values of CQ and HCQ were very similar in pattern, which may be due to structural similarity between CQ and HCQ. Overall, our results indicated that the residue-wise fluctuations of PLpro in apo form were similar to PLpro bound CQ and HCQ as well as PLpro bound reference ligands (PB5 and TTT).

For the entire MDS, H-bond formation is the key indicator of specificity and molecular interactions between the protein and inhibitor complexes. The mean values for H-bonds formed between PLpro and inhibitors were calculated for entire trajectories and are plotted in Figure 6A . Figure  6B . The details of the molecular interactions in all PLpro bound ligand complexes after 10 ns interval are depicted in supplementary Figures S7-S10. As seen in supplementary Figure S10 for CQ complex, H-bonding interaction with Tyr269 with terminal protonated amine and ring Van der Waal interaction with Pro 249, Pro248 and Met209 are conserved. However, in case of HCQ, Tyr269 forms H-bonding with 4-aminoquinoline of CQ-PLpro complex.

Conformational space and transition dynamics of apo PLpro and its inhibitor complexes was further evaluated by the Principle Component Analysis (PCA). The PCA is a statistical calculation that decrease the complexity of MDS trajectory data by extracting only collective motion of Ca backbone atoms whilst preserving most of the significant variations to assess the complex stability. It computes the covariance matrix of positional fluctuations for C-alpha backbone atoms to decipher the dynamics of PLpro-inhibitor interactions. The 2-D projection of the trajectories for two major principal components PC1 and PC3 for SARS-CoV-2 PLpro; SARS-CoV-2 PLpro_PB5; SARS-CoV-2 PLpro_TTT; SARS-CoV-2 PLpro_CQ and SARS-CoV-2 PLpro_HCQ represents different conformations in 2-D space are shown in Figure 7 . The PCA analysis revealed the following observations. First, the 2-D projections of apo PLpro and PLpro-inhibitor were almost identical. Second, the 2-D plot revealed that the SARS-CoV-2 PLpro_TTT ( Figure 7D ) and SARS-CoV-2 PLpro_HCQ ( Figure  7F ) showed higher stability and occupy lesser phase space compared to others. The covariance and 2-D plot analysis also indicated the presence of two well-defined clusters in SARS-CoV-2 PLpro_CQ ( Figure 7E ) and SARS-CoV-2 PLpro_PB5 ( Figure 7C ) whereas SARS-CoV-2 PLpro ( Figure  7B ); SARS-CoV-2 PLpro_TTT ( Figure 7D ) and SARS-CoV-2 PLpro_HCQ ( Figure 7F ) showed only one defined cluster. The positive and negative limits are depicted by the covariance plots where positive values are related to the motion of the atoms occurring along the same direction (correlated), whereas the negative values indicate motion of the atoms in the opposite direction (anti-correlated). Our PCA analysis from MD simulations (50 ns) revealed that the apo SARS-CoV-2 PLpro had more opposite direction (anti-correlated) motion, whilst PLpro-inhibitor complexes had a balance of correlated as well as anti-correlated motion. A plot of eigenvalues calculated from the covariance matrix of backbone fluctuations, plotted in decreasing order versus the respective eigenvector indices for all MD systems is plotted in supplementary Figure S6 . For visualisation of essential dynamics, 50 structures were extracted from each MD simulation projecting the extremely selected eigenvectors. Thus, it was concluded that the presence of inhibitor complex with PLpro indicated more stable and less flexible dynamics.

3.3.5. Binding free energy estimation and energy decomposition of SARS-CoV-2 PLpro complexes with inhibitors

The binding free energy (DG) of PLpro-inhibitor complexes was calculated using the MM-PBSA method for last 30 ns converged frames to give estimated non-bonded interaction energy. A total of 150 frames at every 200 ps from last 30 ns trajectories were taken for the calculation. The estimated value of DG (in kJ mol À1 ) were PB5 (À113.379 ± 1.093), TTT (À75.953 ± 1.133), CQ (À148.178 ± 1.991) and HCQ (À119.678 ± 1.788), respectively (Figure 8a ). The DG value for CQ and HCQ were higher than the reference known crystal structure inhibitor. The values of SARS-CoV-2 PLpro_TTT were lesser than other reference inhibitor PB5. Whilst in the case of CQ and HCQ, the CQ complex with PLpro exhibited highest DG values. Furthermore, the individual component for binding energy, the electrostatic interactions, the Van der Waals, and non-polar solvation energy except the polar solvation energy had favourably contributed to the overall interaction as shown in Figure 8b . For PLpro-inhibitor complexes, the Van der Waal and the non-polar solvation energy were almost identical, whereas the electrostatic energy was higher for HCQ compared to other three complexes.

Inspired by the fact that CQ and HCQ drugs are effective for COVID-19 treatment, we aspired to understand the binding of CQ and HCQ with SARS-CoV-2-PLPro. Binding mode and interaction energies were modelled and comparison was made with reported binders. Based on the results, estimated DG value for PB5, TTT, CQ and HCQ were À113.379 ± 1.093, À75.953 ± 1.133, À148.178 ± 1.991, À119.678 ± 1.788 kJ mol À1 , respectively (Figure 8a ). Of particular note, the DG values for CLQ and HCQ (SARS-CoV-2 PLpro_CQ and SARS-CoV-2 PLpro_HCQ) were higher than the reference known crystal structure inhibitors.

Recently few studies have indicated Spike glycoprotein, RNA-dependent RNA polymerase, Chimeric Receptor-binding domain, Main protease, Non-structural Protein 3, Non-structural Protein 9, ADP-ribose-1 monophosphatase and NTD-N protein as the drug target of CQ and HCQ (Baildya et al., 2020; Mukherjee et al., 2020; Nimgampalle et al., 2020) . However, in present study, we studied PLpro drug target and aimed to characterise the binding mode, binding affinities and key amino acid residues in CQ and HCQ binding to PLpro. Subsequently, they have been compared with known inhibitors and project them for further exploration for designing potent inhibitors in near future. As an outcome of this study, we anticipate that PLPro is the promising target of CQ and HCQ, therefore need to be tested using experimental studies. Moreover, in view of on-going debate around these two molecules, our results will help to promote allied research and drug discovery against the growing epidemic. 

No potential conflict of interest was reported by the authors.

http://orcid.org/0000-0002-1811-2057 Mohd Athar http://orcid.org/0000-0001-6337-1026 P. C Jha http://orcid.org/0000-0002-1709-511X

Repurposing of chloroquine and some clinically approved antiviral drugs as effective therapeutics to prevent cellular entry and replication of coronavirus

GS-5734) is mediated by the viral polymerase and the proofreading exoribonuclease

Docking study of chloroquine and hydroxychloroquine interaction with RNA binding domain of nucleocapsid phospho-protein -An in silico insight into the comparative efficacy of repurposing antiviral drugs

X-ray structural and biological evaluation of a series of potent and highly selective inhibitors of human coronavirus papain-like proteases

The SARS-coronavirus papain-like protease: Structure, function and inhibition by designed antiviral compounds

Inhibitory activity of hydroxychloroquine on COVID-19 main protease: An insight from MD-simulation studies

The papain-like protease of severe acute respiratory syndrome coronavirus has deubiquitinating activity

The 2019-new coronavirus epidemic: Evidence for virus evolution

Molecular dynamics with coupling to an external bath

Dassault Syst emes

Novel 2019 coronavirus structure, mechanism of action, antiviral drug promises and rule out against its treatment

A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: A study of a family cluster

Emerging coronaviruses: Genome structure, replication, and pathogenesis

Role of lopinavir/ritonavir in the treatment of SARS: Initial virological and clinical findings

Chloroquine and hydroxychloroquine as available weapons to fight COVID-19

Structural insights into the interaction of coronavirus papain-like proteases and interferon-stimulated gene product 15 from different species

Particle mesh Ewald: An NÁlog(N) method for Ewald sums in large systems

Peptide folding: When simulation meets experiment

Breakthrough: Chloroquine phosphate has shown apparent efficacy in treatment of COVID-19 associated pneumonia in clinical studies

PH-dependent entry of chikungunya virus into Aedes albopictus cells

The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities

Structurebased design, synthesis, and biological evaluation of a series of novel and reversible inhibitors for the severe acute respiratory syndromecoronavirus papain-like protease

Glide: A new approach for rapid, accurate docking and scoring. 2. Enrichment factors in database screening

LINCS: A linear constraint solver for molecular simulations

The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health -The latest 2019 novel coronavirus outbreak in Wuhan

VMD: Visual molecular dynamics

Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids

A (2014). g_mmpbsa -A GROMACS tool for high-throughput MM-PBSA calculations

Inhibitor recognition specificity of MERS-CoV papain-like protease may differ from that of SARS-CoV

Improved side-chain torsion potentials for the Amber ff99SB protein force field

Structure validation by Ca geometry: u,w and Cb deviation

Protein and ligand preparation: Parameters, protocols, and influence on virtual screening enrichments

Structural insight to hydroxychloroquine-3C-like proteinase complexation from SARS-CoV-2: Inhibitor modelling study through molecular docking and MD-simulation study

Insights into SARS-CoV-2 genome, structure, evolution, pathogenesis and therapies: Structural genomics approach

Screening of chloroquine, hydroxychloroquine and its derivatives for their binding affinity to multiple SARS-CoV-2 protein drug targets

Crystal structure and pair potentials: A molecular-dynamics study

L. donovani XPRT: Molecular characterization and evaluation of inhibitors

Exploring Ruthenium-based organometallic inhibitors against Plasmodium Calcium Dependent Kinase 2 (PfCDPK2): a combined ensemble docking, QM paramterization and molecular dynamics study Funding Sources

Inhibition of amyloid fibril formation of lysozyme by ascorbic acid and a probable mechanism of action

Combined in silico approaches for the identification of novel inhibitors of human islet amyloid polypeptide (hIAPP) fibrillation

UCSF Chimera-a visualization system for exploratory research and analysis

A noncovalent class of papain-like protease/ deubiquitinase inhibitors blocks SARS virus replication

Severe acute respiratory syndrome coronavirus papain-like protease: Structure of a viral deubiquitinating enzyme

Prediction of absolute solvation free energies using molecular dynamics free energy perturbation and the opls force field

ACPYPE -AnteChamber PYthon Parser interfacE

Assessing the performance of MM/PBSA and MM/GBSA methods. 7. Entropy effects on the performance of endpoint binding free energy calculation approaches

GROMACS: Fast, flexible, and free

Chloroquine is a potent inhibitor of SARS coronavirus infection and spread

A novel coronavirus outbreak of global health concern

Calculating protein-ligand binding affinities with MMPBSA: Method and error analysis

Remdesivir and chloroquine effectively inhibit the recently emerged novel coronavirus (2019-nCoV) in vitro

Coronavirus pathogenesis

Improving the physical realism and structural accuracy of protein models by a two-step atomic-level energy minimization

The I-TASSER Suite: Protein structure and function prediction

A systematic review of lopinavir therapy for SARS coronavirus and MERS coronavirus -A possible reference for coronavirus disease-19 treatment option

Discovery of a novel coronavirus associated with the recent pneumonia outbreak in humans and its potential bat origin

A pneumonia outbreak associated with a new coronavirus of probable bat origin

A novel coronavirus from patients with Pneumonia in China