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

The spread of novel coronavirus strain, Severe Acute Respiratory Syndrome 2 (SARS-CoV-2) causes Coronavirus disease (COVID-19) has now spread worldwide and effecting the entire human race. The viral genetic material is transcripted and replicated by 3 C-like protease, as a result, it is an important drug target for COVID-19. Hydroxychloroquine (HCQ) report promising results against this drug target so, we perform molecular docking followed by MD-simulation studies of HCQ and modelled some ligand (Mod-I and Mod-II) molecules with SARS-CoV-2-main protease which reveals the structural organization of the active site residues and presence of a conserve water-mediated catalytic triad that helps in the recognition of Mod-I/II ligand molecules. The study may be helpful to gain a detailed structural insight on the presence of water-mediated catalytic triad which could be useful for inhibitor modelling. Communicated by Ramaswamy H. Sarma.

Introduction

The novel coronavirus outbreak was first reported in Wuhan city of China on December 2019 which causes the Coronavirus disease 2019 (COVID-19) (Sohrabi et al., 2020) is now spread worldwide and effecting the entire human race. The virus strain responsible for this ongoing pandemic is called Severe Acute Respiratory Syndrome 2 (SARS-CoV-2) which is the successor to SARS-CoV-1 and also related to Middle East coronavirus (MERS-CoV) (Petrosillo et al., 2020) having the highest mortality rate among all human coronaviruses (Elfiky & Azzam, 2020). All three viruses have a common host and have originated from bats (Hemida & Ba Abduallah, 2020; Hu et al., 2015). The spread of SARS-CoV-1 (in the year 2002) and MERS-CoV (in the year 2012) to humans is believed to be through civet cats (Wang & Eaton, 2007) and camels (Azhar et al., 2014) unfortunately the mediator for SARS-CoV-2 is still unknown but researchers suspect that ant-eating pangolins may be responsible for this spread as Pangolin-CoV is 91.02% identical to SARS-CoV-2 (T. Zhang et al., 2020).

During this pandemic several drug targets have been identified which includes RNA-dependent RNA polymerase (Anwar et al., 2020; Arya & Dwivedi, 2020; Babadaei, Hasan, Vahdani, et al., 2020; Elfiky, 2020b), Spike glycoprotein (Adeoye et al., 2020; Basit et al., 2020; de Oliveira et al., 2020; Elfiky, 2020a; Sinha et al., 2020; Wahedi et al., 2020), Angiotensin-converting enzyme 2 (ACE2) (Abdelli et al., 2020; Ahmad et al., 2020; Babadaei, Hasan, Bloukh, et al., 2020; Hasan et al., 2020; V. Kumar, Dhanjal, Bhargava, et al., 2020; Veeramachaneni et al., 2020; Wahedi et al., 2020), Nucleocapsid protein (Sarma et al., 2020), Papain-like proteinase (Quimque et al., 2020), 3 C-like proteinase (Das et al., 2020; Enmozhi et al., 2020; Islam et al., 2020; Umesh et al., 2020) and they may be useful for drug development. However, design and search of multi-epitope vaccine candidate molecules have also been initiated in different laboratories (Abraham Peele et al., 2020; Enayatkhani et al., 2020).

The SARS-CoV-2 is a positive-sense single-stranded RNA virus (Hendaus, 2020; Lu et al., 2020) which utilizes non-structural protein (nsp) complexes for the process of transcription and replication in human cells. The replicase polyprotein is processed by viral protease known as Papain like protease (PLP) and 3 C-like protease (3CL protease) or also known as the main protease which is attractive drug targets for the present pandemic COVID-19 (Al-Khafaji et al., 2020; Bhardwaj et al., 2020; Gupta et al., 2020; Gyebi et al., 2020; A. Kumar, Choudhir, Shukla, et al., 2020; D. Kumar, Kumari, Jayaraj, et al., 2020; Mittal et al., 2020; Sk et al., 2020) and is most widely studied target for COVID-19 which also remains conserved among several coronaviruses (CoV), so we select it for our study. The X-ray crystal structures of SARS-CoV-2 3CL-protease were surveyed from RCSB Protein Data Bank (Berman et al., 2000; Burley et al., 2019) which showed the protein consists of three domains. The Domains I (consisting of 1–99 amino acids), Domain II (consisting of 100–183 amino acids), and Domain III (consisting of 200–306 amino acids) in which domain I contributed His41 and domain II contributed Cys145 together forming a catalytic dyad found in entire CoV superfamily. However, the third catalytic residue is missing and the position seems to be occupied by an invariant water molecule (L. Zhang, Lin et al., 2020). The Domain I is connected to Domain III through a connecting loop consisting of 184–199 amino acid residues (Figure 1(a)). The ligand-binding site is comprised of Thr26, His41, Phe140, Asn142, Gly143, Cys145, His163, His164, Met165, Glu166, Gln189 and Thr190 residues as can be seen in 6lu7 crystal structure(Jin et al., 2020) (Figure 1(b)).

Figure 1.

The structural representation of SARS-CoV-2 Mpro (PDB Id:6LU7). a) Domain I, II III, and the connecting loop are shown by yellow, cyan, green and orange colours respectively. This colouring scheme for the protein has been maintained in the entire manuscript. The catalytic dyad (His41 and Cys145) and the invariant water molecule (HOH445) are shown here. b) The active site residues with the N3-inhibitor molecule (carbon atom is shown with white). The carbon atoms of amino acids are shown with their representative domain colour.

The Hydroxychloroquine is an FDA approved drug that has been used to treat malaria(Al-Bari, 2015; Tanenbaum & Tuffanelli, 1980) rheumatoid arthritis(Rynes, 1988), chronic discoid lupus erythematosus(Wahie & Meggitt, 2013), and systemic lupus erythematosus(Chew et al., 2020). During this pandemic emergency, HCQ has been showing promising effect in the treatment of COVID-19 (Boopathi et al., 2020; Hendaus & Jomha, 2020) so, several theoretical and experimental studies are conducted to understand the inhibiting mechanism of this malarial drug to SARS-CoV-2 (Amin & Abbas, 2020; Beura & Chetti, 2020; Beura & Prabhakar, 2020; Devaux et al., 2020; Gautret et al., 2020).

The X-ray crystallographic studies have indicated the inhibitor/drug binding site on 3CL-protease which helps in identification for new potential inhibitor molecules, consequently a large number of theoretical studies are conducted that attempts identification of repurposing drugs for COVID-19 (Aanouz et al., 2020; Choudhury, 2020; Elmezayen, Al-Obaidi, Şahin, & Yelekçi, 2020; Joshi et al., 2020; M. T. Khan et al., 2020; R. J. Khan et al., 2020; S. A. Khan, Zia, Ashraf, Uddin, & Ul-Haq, 2020; V. Kumar, Dhanjal, Kaul, Wadhwa, & Sundar, 2020; Lobo-Galo, Terrazas-López, Martínez-Martínez, & Díaz-Sánchez, 2020; Mahanta et al., 2020; Muralidharan, Sakthivel, Velmurugan, & Gromiha, 2020; Pant, Singh, Ravichandiran, Murty, & Srivastava, 2020). In our study we carefully investigated the catalytic and inhibitor binding site of 3CL-protease in the crystal and MD-simulated structures. We also report the binding mechanism and stabilization of Hydroxychloroquine with SARS-CoV-2 main protease through non-covalent interactions as well as attempt has been made to model some new structural analogs (Mod I and I/, Mod II and II/) and their recognition dynamics to conserved water-mediated catalytic triad.

Material and methods

The X-ray structures of some 3CL-protease from SARS-CoV-2 have been surveyed which were selected from RCSB Protein Data Bank (Supplementry material Table S1). The position and hydrogen bonding interaction of catalytic residues were examined using the UCSF Chimera program(Pettersen et al., 2004).

Table 1.

The MM-PBSA binding free energy calculation for HCQ, Mod I, and Mod II-MPro complex structures during MD-simulation. All the given energies are in kcal/mol.

Complex	Time (ns)	ΔEElectrostatic	ΔEVan der Waal	ΔEpolar	SASA	ΔEbinding
HCQ-MPro	20	−130.55 +/- 3.47	−35.81 +/- 3.74	36.24 +/- 2.95	−40.70 +/- 3.65	−170.82 +/- 3.45
	40	−137.29 +/- 1.38	−34.66 +/- 1.33	33.57 +/- 1.16	−39.44 +/- 1.30	−177.82 +/- 1.29
	60	−151.44 +/- 3.22	−32.61 +/- 3.12	29.57 +/- 2.83	−37.27 +/- 3.10	−191.75 +/- 3.06
	80	−150.93 +/- 2.72	−33.08 +/- 2.87	31.86 +/- 1.89	−37.74 +/- 2.84	−189.89 +/- 2.58
	100	−153.59 +/- 4.95	−32.72 +/- 2.59	30.96 +/- 1.65	−37.34 +/- 2.57	−192.69 +/- 2.94
Mod I-MPro	20	−115.93 +/- 3.93	−36.14 +/- 3.98	39.10 +/- 0.30	−40.84 +/- 4.06	−153.81 +/- 3.06
	40	−124.83 +/- 2.16	−39.37 +/- 1.99	42.01 +/- 1.54	−44.48 +/- 1.96	−166.67 +/- 1.91
	60	−125.62 +/- 2.54	−39.60 +/- 4.23	46.33 +/- 2.55	−44.70 +/- 4.21	−163.59 +/- 3.38
	80	−125.24 +/- 4.47	−39.10 +/- 2.34	45.54 +/- 3.30	−44.22 +/- 2.34	−163.02 +/- 3.11
	100	−135.22 +/-2.66	−35.09 +/- 2.09	39.68 +/- 2.44	−40.04 +/-2.06	−170.67 +/- 2.31
Mod II-MPro	20	−148.43 +/- 6.29	−48.22 +/- 3.38	51.47 +/- 4.73	−54.02 +/- 3.42	−199.2 +/- 4.45
	40	−157.48 +/- 2.42	−45.61 +/- 1.89	43.25 +/- 0.98	−51.15 +/- 1.89	−210.99 +/- 1.79
	60	−157.20 +/- 1.06	−45.32 +/- 2.80	43.34 +/- 0.91	−50.92 +/- 2.76	−210.1 +/- 1.88
	80	−155.53 +/- 4.03	−46.20 +/- 2.80	42.28 +/- 0.85	−51.86 +/- 2.79	−211.31 +/- 2.61
	100	−157.45 +/- 1.00	−45.88 +/- 3.41	41.87 +/- 2.64	−51.63 +/- 3.38	−213.09 +/- 2.60

Ligand preparation

The ligand (HCQ, Mod I, Mod I/, Mod II, Mod II/) structures (Figure 2) were built by using the Gaussview program(Dennington et al., 2009). The structures were then subjected to geometry optimization with B3LYP/6-31G level of theory, solvation energies were added using the conductor-like Polarizable Continuum Models (CPCM). The partial atomic charges of the atoms were retained during classical MD-simulation. All the QM calculations were performed by using the Gaussian 09 package(Frisch et al. (2009) Gaussian 09 A.02).

Figure 2.

Chemical structures of HCQ and modified ligand molecules included in this study.

Protein–ligand docking

The 3CL-protease has one ligand binding site which has been indicated in the X-ray crystallographic studies. Ligand–receptor docking was performed using Autodock Vina v.1.1.1 (Trott & Olson, 2009). The 6LU7 X-ray structure (excluding water and ligand molecules) was considered as a receptor. Two PDBQT files were generated for the receptor protein, one for rigid portion and another for flexible (His41 and Cys145) side chains using AutoDock Tools v.1.5.4 (Morris et al., 2009) by assigning Kollman united atom charges (Weiner et al., 1984). The QM optimized structures of ligand molecules (HCQ and Mod-I, I/, II, II/) were converted into PDBQT files after including their partial atomic charges using the Gasteiger method (Gasteiger & Marsili, 1980). Grid point spacing was set at 1 Å and 20 grid points were taken in each direction. As the location of the ligand binding site was already known (PDB id: 6LU7), so grid box was centered at that site. Vina automatically calculated the grid map for searching. All other docking parameters were assigned to their default values. The docked structures were visually inspected and selected for further work. The ligand–receptor docking was also validated with DockThor (A receptor-ligand docking program) server (De Magalhães et al., 2014; Santos et al., 2020).

Molecular dynamics (MD) simulation

Molecular dynamics simulation was performed using NAMD v.2.12 (Kalé et al., 1999; Phillips et al., 2005) with the CHARMM36 force field(Brooks et al., 1983; Huang & MacKerell, 2013; MacKerell et al., 1998) . Then each structure was converted to Protein Structure File (PSF) by Automatic PSF Generation Plug-in within VMD program v. 1.9.2 (Humphrey et al., 1996). Subsequent energy minimization was performed by the conjugate gradient method. The process was conducted in two successive stages; initial energy minimization was performed for 2000 steps by fixing the backbone atoms, followed by a final minimization for 5000 steps that were carried out for all atoms of the system to ensure the removal of any residual steric clashes. Then the energy minimized structures were simulated at constant temperature (310 K) and pressure (1 atm) by Langevin dynamics (Gullingsrud et al., 2001) using periodic boundary condition. The Particle Mesh Ewald method was applied for full-electrostatics and the Nose–Hoover Langevin piston method was used to control the pressure and dynamical properties of the barostat(Feller et al., 1995). All-atom molecular dynamics simulation for 100 ns was carried out for HCQ, Mod I, and Mod II docked 3CL-protease structure. The atomic coordinates of simulated structures were recorded at every 2 ps for further analysis.

MM-PBSA calculation

The ligand-receptor binding energies were calculated by Molecular Mechanics-Poisson Boltzmann Surface Area (MM-PBSA) method using Calculation of Free Energy (CaFE) program(Liu & Hou, 2016). It is capable of handling the simulation trajectories generated by various force field and is powered by VMD and NAMD. Three different energy terms are calculated by the MM-PBSA method which includes energy difference between the complex and separate receptor and ligand molecule in the gas phase from NAMD. Then, the polar solvation free energy was calculated using the APBS program(Baker et al., 2001) and finally, the difference of solvent-accessible surface area (SASA) is measured and the nonpolar solvation free energy was estimated. The binding energies are calculated by using the following equation:

The binding energies of the ligand-receptor complexes were calculated at 20 ns interval over the simulation trajectories.

Results and discussion

The analysis of 3CL-protease crystal structures from SARS-CoV-2 showed that the His 41 of Domain I lies close to Cys145 of Domain II which together form the catalytic dyad. An invariant water molecule, HOH445 (W1) is positioned within hydrogen bond distance to His41ND1, Arg40NB, His164ND1, and Asp187OD2 atoms (Figure 3).

Figure 3.

The interaction of His41ND1, Arg40NB, His164ND1, and Asp187OD2 residues with an invariant water molecule through hydrogen bond interaction. The hydrogen bonds are shown by dotted lines.

The Root Mean Square Deviation (RMSD) of the protein backbone (Cα) atoms (Figure 4) show that the Mod I/Mod II complexed MPro structures are more stable than the HCQ complexed structure whereas, the Root Mean Square Fluctuation (RMSF) curve shows the active site residues in all the three complexes are stable with small fluctuations at N and C-terminal region of the protein. All the ligand-protein complexes show similar deviation and fluctuation validating the consistency of the simulation trajectories.

Figure 4.

The Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) curves of the protein backbone (Cα) atoms during MD-simulation.

Molecular docking and simulation studies of hydroxychloroquine (HCQ) with 3CLpro

The optimized structure of HCQ was docked in 6lu7 crystal structure, the ligand molecule binds to a groove on the surface of the protein with binding energy ∼-5.4 kcal/mol. The catalytic dyad (His41 and Cys145) interacts with the si-face (under the aromatic ring) of the ligand molecule through non-covalent interaction. The His41 endowed parallel displaced π‧‧‧π interaction with aromatic ring 1 of HCQ and the Cys145 shown S-H‧‧‧π interaction with aromatic ring 2 of HCQ. The Thr24OG1 located at the entrance of the inhibitor binding site is stabilized by the O21 hydroxyl group of HCQ through hydrogen bond interaction (Supplementry material Figure S1).

During simulation, the hydrogen bond interaction between Cys145SG and His41NE2 is ranging from 3.1 to 3.3 Å and the average H-bond interaction of W1 water center to His41ND1, Arg40NB, His164ND1, and Asp187OD2 atoms are observed to be ∼ 3.0–3.3, 3.1–3.4, 2.8–3.4, 2.7–3.2 Å respectively (Figure 5).

Figure 5.

Hydrogen bond interaction between different residues in HCQ-SARS-CoV-2 Mpro during MD-simulation.

A distorted trigonal pyramidal geometry has been observed around the W1 water center where Arg40NB and His164 at the apex position and His41, His164, and Asp187 are present at the three corners of the plane (Figure 6) with ∼100% residential frequency. In this HCQ complex structure, the Cys145SG atom is observed to lie near C8 and C9 carbon centers with ∼3.88 and 3.72 Å distance respectively suggesting these two centers may be important for inhibitor modeling so, these carbon centers were modified by substituting with different functional groups for easy nucleophilic attack by Cys145 catalytic residue.

Figure 6.

The hydrogen bonding geometry around W1 water center involving His41, Arg40, His164, and Asp187 residues during MD-simulation.

Molecular docking and simulation studies of modified inhibitor molecules with 3CL-protease

Observing the structural arrangement and interaction of His41 and Cys145 with HCQ, some ligand molecules were modelled (Figure 2) and docked into the inhibitor binding cavity of the 6lu7 crystal structure. The docked structures of modified ligand molecules (Mod I, I/, II, and II/) win the 6lu7 crystal structure are shown in Supplementry material Figure S2.

The binding energies (B.E.) for Mod I and Mod II are observed to be −6.1 and −5.5 kcal/mol whereas that of Mod I/and Mod II/are −5.4 and −5.3 kcal/mol respectively. The results indicate that the modification of HCQ at the C8 position is energetically favourable compared to the C9-position. Moreover, we may say the formation of carbocation at C8-position of Modified ligand molecules would be more stable than at C9-position due to –I effect of N7-atom in the aromatic ring of the ligand molecule so, Mod I and Mod II structures were selected for further MD-simulation studies.

The molecular dynamics simulation of Mod I and Mod II complexed 6lu7 structure shows that the Cys145SG lies in close proximity to O25 atom of the ligand with average distance ranging from 3.2 to 3.7 Å for Mod I and 3.5 to 4.0 Å for Mod II complexed structure (Figure 7 and 8). The hydrogen bond interaction between the catalytic dyad Cys145SG‧‧‧His41NE2 varies from 3.2 to 3.5 Å. The W1 water molecule forms strong hydrogen bond interaction with His41ND1 and Asp187OD2 thus bridging the gap between these two residues. The His41ND1‧‧‧W1 distance ranges from ∼3.1 to 3.4 Å and 3.0 to 3.3 Å whereas that of W1‧‧‧Asp187OD2 ranges from 2.7 to 2.9 and 2.7 to 3.0 Å for the respective Mod I and Mod II complex structures (Figure 7 and 8).

Figure 7.

Hydrogen bonding interaction between the different residues in Mod-I-SARS-CoV-2 Mpro during MD-simulation.

Figure 8.

Hydrogen bond interaction between different residues in Mod-II-SARS-CoV-2 Mpro during MD-simulation.

In classical proteases (like chymotrypsin, trypsin or elastase), the catalytic triad is made up of an acid, base and a nucleophile but in coronavirus main proteases (SARS-CoV-2 Mpro (PDB Id 6lu7), SARS-CoV Mpro (PDB Id 2zu5) and MERS-CoV Mpro (PDB Id 5wkk)) presence of water-mediated Aspartic acid residue along with base (His) and nucleophile (Cys) seem to be unique and reflecting the formation of water-mediated catalytic triad so, the occupation of W1 water center adjacent to the base (His41) and its recognition to Asp187 through Asp187OD2‧‧‧W1‧‧‧His41NE2 interaction could be an evolutionary change in the catalytic site of CoV-Mpro for unknown reasons.

During simulation, the recognition of Asp187 to Mod I/II is observed to made through Asp187OD2‧‧‧W1‧‧‧His41ND1/NE2‧‧‧Cys145SG‧‧‧Mod I/II (O25) path. The carbonyl oxygen (O25) is stabilized by amide backbone atoms of Cys145, Ser144, and Gly143 residues (Figure 9).

Figure 9.

Recognition of Mod-I/II to Cys145, His41, Asp187 residues and W1(HOH445) water center.

MM-PBSA binding free energy

The binding free energies of the ligand-receptor complexes by the MM-PBSA method (Table 1) are calculated. The final binding energy is the cumulative sum of Electrostatic, Van der Waal, Polar solvation, and SASA energies which show that the Mod II structure has better binding energies when compared to Mod I and HCQ molecules. The Mod I structure shows low binding and electrostatic energies compared to the HCQ complex structure. The Mod II complexed structure shows high binding energy of ∼−213.09 kcal/mol at 80–100ns time interval.

So, analyses of crystal and ligand-bound simulated structures from SARS-CoV-2 Mpro a mechanism has been proposed which may shed some light on the mechanistic role of the above residues in ligand complexed Mpro structures (Figure 10).

Figure 10.

Proposed mechanism of action in SARS-CoV-2 Mpro.

In SARS-CoV-2 Mpro the catalytic triad is made up of Cys145, His41, and Asp187 residues, a catalytic water (HOH445) which is associated with His41 and Asp187 through hydrogen bond interaction. Initially, His41 abstracts a proton from Cys145 making it a better nucleophile (cysteine thiyl radical, CysS-) for attack on the carbonyl carbon (C24) of the ligand molecule. The Asp187 and HOH445 anchor the His41 molecule to its correct conformation during the formation of transition state complexes. The reaction is followed by the formation of a tetrahedral intermediate where the oxyanion group of the intermediate is stabilized by amide backbone atoms of Gly143, Cys145, and Ser144 residues. Finally, the ROH product molecule is released from the reaction medium by the formation of a stable Mod-lig-Cys145 adduct.

Conclusion

The Molecular docking followed by MD-simulation of the HCQ complexed structure of SARS-CoV-2 Mpro reveals the binding of this ligand molecule on the surface of the protein structure near to the catalytic dyad (His41 and Cys145). The Cys145 is observed to lie near the C8 and C9 centers of the ligand molecule which encouraged us to modify these centers by addition of –COOH and –COOCH3 functional groups. The modified ligand molecules (Mod I and Mod II) are found to have higher binding energy towards the inhibitor binding domain compared to their counter partners (Mod I/and Mod II/). During the simulation of Modified ligand molecules with SARS-CoV-2 Mpro, the Cys145 is observed to be perfectly positioned near the C24 atom of ligand molecules thus suggesting the possibility of a nucleophilic attack to this site during catalysis. The W1 water center hydrogen-bonded to His41ND1 seems to be crucial for the recognition of ligand molecules to catalytic dyad and Asp187 (located on the connecting loop between Domain I and Domain III) through Mod-I/II(O25)‧‧‧Cys145‧‧‧His41NE2/ND1‧‧‧W1‧‧‧Asp187OD2 interaction. So, from the results, it may be presumed that the Asp187 may function as the third catalytic residue where W1 water center maintains the structural recognition of this catalytic triad.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.