Computational investigation of binding of chloroquinone and hydroxychloroquinone against PLPro of SARS-CoV-2.

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.

Introduction

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; Wang 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; Zhou et al., 2020).

The causative agent for SARS-CoV-2, coronavirus belongs to subgenus Sarbecovirus, Coronaviridae family and Orthocoronavirinae subfamily. It is a β‐coronavirus of 2B (β-coronavirus), however the group is non-zoonotic till date unlike earlier coronaviruses (Hui et al., 2020; 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 (Zhou et al., 2020; 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 (Zhou et al., 2020). 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) (Chen et al., 2020). 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 (Zhu et al., 2020). 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áez-Santos et al., 2014, 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; Wang 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; Wang et al., 2020) and Favipiravir as inhibitors of RNA-dependent RNA polymerase (Wang et al., 2020), 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.

Materials and methods

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).

Ligand preparation and molecular docking

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 multi-conformation 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.

Molecular dynamics (MD) simulations of SARS-CoV-2 PLpro and inhibitors

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, 2018). 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). After the energy minimisation, equilibration with position restraint was carried out under NVT (constant number [N], constant volume [V] and constant temperature [T]) and NPT (constant number [N], constant pressure [P] and constant temperature [T]) ensemble for 1 ns each. For NVT equilibration, modified Berendsen thermostat algorithm (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) (Patel et al., 2018). For molecular visualization and other analysis Discovery studio was utilized (BIOVIA 2020).

Clustering of conformations and PCA/essential dynamics

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 Cα 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 Cα 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.

Estimation of binding free energy using MM/PBSA

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.

Results and discussion

Model building and validation of SARS-CoV-2 PLpro

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.

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.

Figure 2.

Three-dimensional structure of SARS-CoV-2 PLpro. The SARS-CoV-2 PLpro model is shown in ribbon model coloured by its secondary structure masked with surface rendering. Domain assignment is marked in the model. Amino acid sites where differ between SARS-CoV-2 and SARS-CoV are coloured in blue. The catalytic triad residues are labelled in magenta, whilst residues forming the zinc-finger motif are labelled with orange colour.

Molecular docking of CQ and HCQ

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 H-bonds, one with N of Tyr269 and another with protonated NH with Asp165. TTT also interacted with Asp165 through H-bonding 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 H-bond 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.

Table 1.

Docking scores and list of major interactions.

Ligand	G score	Interactions
PB85 (4ovz)	−6.97	Hydrophobic: Pro249, Pro248, Tyr269 (pi-pi), Met209, Gln270, Leu263, Tyr274H-bond: Asp165, Tyr265Ionisable: Gln270, Asn260
TTT(3e9s)	−5.37	Hydrophobic: Gly164, Leu163, Tyr274, Met209, PRO249, Pro248, Tyr265, tyr269(pi-pi)H-bond: Glu168, Asp165Ionisable: Asn268, Arg167, Thr302
CLQ	−4.72	Hydrophobic: Met209, Pro249, Pro248, Thr302, Tyr274, Asp165, Gly164, Leu163, Tyr265, Tyr269 (pi-pi)H-bond: Tyr269Ionisable: Asn268, Gln270, Asp165
HQ_mode1	−4.09	Hydrophobic: Tyr274, Tyr269, Pro249, Pro248, Ala247, Val166, arg167 (cation-pi), Met209H-bond: Asp165, Tyr265, Arg167Ionisable: Glu168, Thr302, Asp303, Ser246, Arg167, Glu168
HQ_mode2	–	Hydrophobic: Tyr265, Pro248, Tyr274, Ley163, Gly164, Tyr269(pi-pi)Ionisable: Asp165, Asn268, Gln270, Thr302H-bond: Tyr269

Molecular dynamics simulations of SARS-CoV-2 PLpro and inhibitors

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.

Binding mode of HCQ (yellow), TTT (green), CQ (pink) and PB5 (cyan) ligands in complex with PLpro.

Stability of SARS-CoV-2 PLpro complexes with inhibitors

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 Cα backbone and ligand for the entire 50 ns simulations. The RMSD measures the structural variation between the Cα 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.

Figure 4.

Molecular Dynamics Simulations of SARS-CoV-2 PLpro in apo and inhibitor bound complexes, computing the deviation (nm) versus function of time (50 ns): (A) RMSD of the protein Cα backbone atoms of SARS-CoV-2 PLpro (black); SARS-CoV-2 PLpro_PB5 (red); SARS-CoV-2 PLpro_TTT (green); SARS-CoV-2 PLpro_CQ (blue) and SARS-CoV-2 PLpro_HCQ (yellow). (B) RMSD of the inhibitor atoms of docked complexes where colour representation is PB5 (black); TTT (red); CQ (green) and HCQ (blue).

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 Cα 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.

Residue-wise fluctuations of SARS-CoV-2 PLpro with inhibitors

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).

Figure 5.

Residue-wise RMSF deviations (nm) of SARS-CoV-2 PLpro in apo and inhibitor bound complexes: (A) RMSF deviation plot of the protein Cα backbone atoms of SARS-CoV-2 PLpro (black); SARS-CoV-2 PLpro_PB5 (red); SARS-CoV-2 PLpro_TTT (green); SARS-CoV-2 PLpro_CQ (blue) and SARS-CoV-2 PLpro_HCQ (yellow). (B) RMSF deviation plot of the inhibitor atoms of docked complexes where colour representation is PB5 (black); TTT (red); CQ (green) and HCQ (blue).

Molecular interactions of SARS-CoV-2 PLpro with inhibitors

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. The average H-bonds formed for PB5; TTT; CQ and HCQ were 1.88 ± 0.32, 1.05 ± 0.65, 0.98 ± 0.11, and 2.95 ± 0.70, respectively. Owing to the additional hydroxyl group, HCQ complex showed the highest number of H-bond formation, whilst the TTT and CQ had a similar number of H-bonds throughout the trajectories. The intra H-bonds in PLpro were also computed for entire trajectories for apo as well as PLpro-inhibitor complexes. The average intra H-bonds formed 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 were 229.19 ± 7.45, 232.1 ± 6.84, 232.24 ± 6.95, 231.84 ± 6.56, and 232.22 ± 6.87, respectively. The values for apo PLpro, as well as PLpro-inhibitor complexes, are almost identical, indicating a stable MDS system in the absence and presence of inhibitors 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.

Figure 6.

Inter and Intra hydrogen bonds of SARS-CoV-2 PLpro in apo and inhibitor bound complexes: (A) Number of hydrogen bonds computed versus function of time (50 ns) between protein and inhibitor in SARS-CoV-2 PLpro_PB5 (black); SARS-CoV-2 PLpro_TTT (red); SARS-CoV-2 PLpro_CQ (green) and SARS-CoV-2 PLpro_HCQ (blue). (B) Intra H-bond formation plot within protein of SARS-CoV-2 PLpro (black); SARS-CoV-2 PLpro_PB5 (red); SARS-CoV-2 PLpro_TTT (green); SARS-CoV-2 PLpro_CQ (blue) and SARS-CoV-2 PLpro_HCQ (yellow).

Sars-CoV-2 PLpro conformation clustering and principle component analysis (PCA)

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 Cα 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.

Figure 7.

PCA 2-D projection scatter plot of SARS-CoV-2 PLpro in apo and inhibitor bound complexes: (A) Overlay of 2-D scatter plot projection the motion of the proteins in phase space for the two principle components, PC1 and PC3 derived from all MD simulation setup. Panel B, C, D, E and F represents individual 2-D plots of 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 respectively. For all the panels colour representation is SARS-CoV-2 PLpro (black), SARS-CoV-2 PLpro_PB5 (red), SARS-CoV-2 PLpro_TTT (green), SARS-CoV-2 PLpro_CQ (blue) and SARS-CoV-2 PLpro_HCQ (yellow).

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

The binding free energy (ΔG) 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 ΔG (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 ΔG 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 ΔG 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.

Figure 8.

MM-PBSA Calculation for binding free energy. (A) The total binding free energy for all the SARS-CoV-2 PLpro-inhibitor complexes calculated for last 30 ns stable trajectory for a total of 150 frames, each at 200 ps interval. (B) Representative contributions of each energy component for binding free energy for all the SARS-CoV-2 PLpro with respective inhibitor complexes.

Conclusions

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 ΔG 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 ΔG 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.

Click here for additional data file.