Evaluation of green tea polyphenols as novel corona virus (SARS CoV-2) main protease (Mpro) inhibitors - an in silico docking and molecular dynamics simulation study.

Coronavirus disease 2019 (COVID-19) is a viral respiratory disease which caused global health emergency and announced as pandemic disease by World Health Organization. Lack of specific drug molecules or treatment strategy against this disease makes it more devastating. Thus, there is an urgent need of effective drug molecules to fight against COVID-19. The main protease (Mpro) of SARS CoV-2, a key component of this viral replication, is considered as a prime target for anti-COVID-19 drug development. In order to find potent Mpro inhibitors, we have selected eight polyphenols from green tea, as these are already known to exert antiviral activity against many RNA viruses. We have elucidated the binding affinities and binding modes between these polyphenols including a well-known Mpro inhibitor N3 (having binding affinity -7.0 kcal/mol) and Mpro using molecular docking studies. All eight polyphenols exhibit good binding affinity toward Mpro (-7.1 to -9.0 kcal/mol). However, only three polyphenols (epigallocatechin gallate, epicatechingallate and gallocatechin-3-gallate) interact strongly with one or both catalytic residues (His41 and Cys145) of Mpro. Molecular dynamics simulations (100 ns) on these three Mpro-polyphenol systems further reveal that these complexes are highly stable, experience less conformational fluctuations and share similar degree of compactness. Estimation of total number of intermolecular H-bond and MM-GBSA analysis affirm the stability of these three Mpro-polyphenol complexes. Pharmacokinetic analysis additionally suggested that these polyphenols possess favorable drug-likeness characteristics. Altogether, our study shows that these three polyphenols can be used as potential inhibitors against SARS CoV-2 Mpro and are promising drug candidates for COVID-19 treatment.

Introduction

The novel corona virus disease (COVID-19) is spreading in the whole world causing deaths of 308,962 people till the 16th of May 2020 (https://www.worldometers.info/coronavirus/). The World Health Organization (WHO) has already declared COVID-19 as a pandemic disease (Cucinotta & Vanelli, 2020). Symptoms of this disease mainly include fever, cough, sore throat, runny nose and difficulty in breathing (Chen et al., 2020; Ren et al., 2020; Zhu et al., 2020). The disease was initially originated from Wuhan, China and expanded to nearly 215 countries (Zhu et al., 2020). USA, Spain and UK are among the most affected countries due to COVID-19 disease. In the USA alone, a total number of 1,484,287 people are infected and 88,507 people died due to COVID-19. In Spain and UK, 274,367 and 236,711 people are infected and 27,459 and 33,998 people have died there, respectively. This viral disease has not only challenged human health but also enormously affected the global economy.

The etiological agent of COVID-19 is severe acute respiratory syndrome corona virus-2 (SARS CoV-2) which belongs to the genus β-coronavirus (Zheng, 2020). Before SARS CoV-2, two more epidemic diseases were caused by corona viruses namely severe acute respiratory syndrome (SARS) and middle east respiratory syndrome (MERS) (Drosten et al., 2003; Zaki et al., 2012). In 2003, SARS had killed ∼750 people, whereas MERS was the reason behind the death of ∼860 people in 2012 (Mahase, 2020). On the contrary, the mortality rate due to SARS CoV-2 is exceptionally high (Mahase, 2020). SARS CoV-2 is a non-segmented enveloped positive-sense single-stranded RNA virus with ∼29.9 kb genome size (Wu et al., 2020; Zhu et al., 2020). The virus is thought to be originated from the bat and transmitted to humans by other sources (Zheng, 2020). Later on, human to human transmission of this disease is confirmed (Chan et al., 2020). Till date, no drug molecules or specific therapies have been developed to combat COVID-19. Considering the risk factors associated with this viral infection, an effective drug molecule is urgently required for effective treatment and to limit the transmission of this disease.

Many efforts involve the targeting of nonstructural protein, spike protein, the RNA-dependent RNA polymerase (RdRp) and the angiotensin-converting enzyme II (ACE2) entry receptor for the anti-COVID-19 drug development (Abdelli et al., 2020; Babadaei et al., 2020; Basit et al., 2020; Elfiky, 2020a, 2020b; Hasan et al., 2020; Nejadi Babadaei et al., 2020; Sinha et al., 2020; Wahedi et al., 2020). Targeting the components responsible for replication of SARS CoV-2 may be also a good strategy to identify effective drugs for the treatment of COVID-19. Like SARS and MERS corona viruses, the SARS CoV-2 genome also contains two open reading frames namely ORF1a and ORF1ab (Boopathi et al., 2020). These two ORFs help in translating two overlapping viral polyproteins pp1a and pp1ab required for viral replication and transcription (Grum-Tokars et al., 2008; Marra et al., 2003; Thiel et al., 2003). The functional polypeptides are released from the pp1a and pp1ab polyproteins by proteolytic processing by the papain-like proteinase (PLpro) and the 3C-like protease (3CLpro) (Grum-Tokars et al., 2008; Marra et al., 2003; Thiel et al., 2003). PLpro and 3CLpro are responsible for the cleaving of three sites and 11 sites, respectively within the viral genome (Harcourt et al., 2004; Thiel et al., 2003). As 3CLpro cleaves most of the sites of the polypeptide, it is also known as the main protease or Mpro. Different studies have identified Mpro as a cysteine protease with a Cys-His catalytic dyad (His41 and Cys145) in the active site of the protease (Blanchard et al., 2004; Dai et al., 2020; Jin et al., 2020; Osman et al., 2020). As Mpro plays a vital role in polyprotein processing and virus maturation, it is considered to be an important target for designing antiviral drugs against SARS CoV-2 (Anand et al., 2003; Yan et al., 2003). In addition, there is no human homolog of Mpro which makes it an ideal antiviral target (Kim et al., 2016). Apart from this, the high-resolution crystal structure of Mpro protease along with its inhibitor is recently made available, which has eased the way in designing structure-based Mpro-specific inhibitors to combat this viral disease (Jin et al., 2020; Yang et al., 2003). Utilizing computational approach, investigators have identified many small molecules including HIV and malaria drugs as SARS CoV-2 protease inhibitor (Adeoye et al., 2020; Khan et al., 2020; Muralidharan et al., 2020). Interestingly, many antiviral phytochemicals, bioactive compounds from Moroccan medicinal plants (Crocin, Digitoxigenin and β-Eudesmol) and chemical compounds from Indian spices (Carnosol, Rosmanol and Arjunglucoside-I) are proposed to be effective SARS CoV-2 Mpro inhibitors (Aanouz et al., 2020; Das et al., 2020; Enmozhi et al., 2020; Gyebi et al., 2020; Islam et al., 2020; Joshi et al., 2020; Umesh et al., 2020).

Plant-derived natural polyphenols are well known in preventing a wide range of diseases which include viral diseases as well. Recently, Purohit and coworkers have revealed that many bioactive molecules including a polymerized polyphenol (Oolonghomobisflavan-A) from tea plant (Camellia sinensis L.) act as effective SARS CoV-2 Mpro inhibitors (Bhardwaj et al., 2020). Green tea (C. sinensis) also contains eight native monomeric catechins or polyphenolic compounds and those are epigallocatechin gallate (EGCG), epigallocatechin (EGC), epicatechin gallate (ECG), epicatechin (EC), gallocatechin-3-gallate (GCG), gallocatechin (GC), catechin gallate (CG) and catechin (C) (Ai et al., 2019). These green tea polyphenols/catechins are promising compounds in exhibiting antiviral activities. They show antiviral activity against a wide range of human viruses including influenza, hepatitis B, hepatitis C, herpes simplex virus and HIV (Calland et al., 2012; Fassina et al., 2002; Ide et al., 2016; Lyu et al., 2005; Xu et al., 2008). These polyphenols are even active against dengue virus (DENV), Chikungunya virus (CHIKV) and Zika virus (ZIKV) (Carneiro et al., 2016; Ismail & Jusoh, 2017; Mahajan et al., 2020; Weber et al., 2015). But whether these polyphenols show any antiviral activity against SARS CoV-2 by inhibiting the enzymatic activity of its Mpro is far from clear. In this current study, we have undertaken a thorough attempt to find out the potentiality of these eight green tea polyphenols/catechins to inhibit the Mpro using in silico docking studies, molecular dynamics (MD) simulations and binding free energy calculations. This study provides three green tea polyphenols (EGCG, ECG and GCG) which may be served as inhibitors against SARS CoV-2 Mpro.

Materials and methods

Preparation of the ligands

The structures of eight green tea catechins or polyphenols [EGCG, EGC, ECG, EC, GCG, GC, CG and C] were downloaded from PubChem crystal database server in SDF format (https://pubchem.ncbi.nlm.nih.gov). Then all the sdf files were converted to pdb files using pymol software and each of the polyphenol structures were optimized with B3LYP/6-31G* type of basis set by using Gaussian09 software (Frisch & Clemente). Following optimization, each of these was loaded to AutoDock Tools for the final preparation of ligands. Each of the optimized polyphenol structure was then inserted to AutoDock Tools and standard processes were used to obtain the pdbqt files.

Preparation of Mpro

The crystal structure of the SARS CoV-2 Mpro was taken from the RCSB Protein Data Bank (http://www.rcsb.org) (PDB ID: 6LU7) (Jin et al., 2020). Then we checked whether there were any improper bonds, missing hydrogens, side chain anomalies present or not. After correcting all those aspects the structure file was inserted into AutoDock Tools and standard procedures were followed to get pdbqt file format of Mpro (Morris et al., 2008, 2009).

Molecular docking

The docking of Mpro with the eight polyphenols of green tea was performed with the aid of AutoDock Vina, which is now widely used as a successor of AutoDock Tools (Morris et al., 2008, 2009). The binding affinities of polyphenols–Mpro were determined and analyzed using the same software. By seeing the position of active site region, the center of the grid box was chosen to be at X: −11.75, Y: 15.135, Z: 68.856 with a suitable gridbox volume where the ligands can easily be fitted and which covers the entire active site pocket. The available crystal structure (PDB: 6LU7) with a bound N3 molecule was used as a reference to calibrate and optimize the docking procedures. Following calibration and optimization, the same grid box size and other parameters were used for docking studies of all the eight polyphenols and the entire set up was run to obtain different docked conformations. Amongst them, the best suited conformations with lowest root mean square deviation (RMSD) values were selected to calculate the binding energetic between Mpro and green tea polyphenols. The output from AutoDock Vina was rendered with PyMOL (DeLano, 2002).

MD simulation

All the MD simulations were carried out by the GROningen MAchine for Chemical Simulations (GROMACS) version 5.1.2 (Abraham et al., 2015). The OPLS-AA/L force field and TIP3P water model embedded in GROMACS were used for all the MD simulations (Kaminski et al., 2001). The topology file of Mpro protein was prepared by the GROMACS, while the ligand topologies were obtained from the LigParGen (a server from the Jorgensen group to produce OPLS topologies). In order to satisfy minimum image conventions the system was initially accommodated in a cubic box with a distance of 1.4 nm between protein complex and the box. All bond lengths of protein and polyphenols were constrained using the LINCS algorithm while water molecules were restrained by SETTLE algorithm (Hess et al., 1997; Miyamoto & Kollman, 1992). A total of 30105, 30065, 30073 and 30073 water molecules were added to a cubic simulation box containing the unligated Mpro, Mpro–ECG, Mpro–EGCG, Mpro–GCG complexes, respectively. Each system was energy-minimized using steepest descent algorithm and equilibrated to achieve the appropriate volume. Short- and long-range non-bonded interactions were calculated by applying twin-range cutoffs of 0.9 and 1.4 nm, respectively. The leap frog algorithm with time step 2 fs was used for integrating the equation of motion and the neighbor list was updated at every five steps. Long-range electrostatics is treated using the Particle Mesh Ewald method, with a Fourier grid spacing of 0.16 (Darden et al., 1993). Periodic boundary conditions were applied in all three directions. Equilibration of the systems was carried out in two main stages. First, the system was allowed to heat gradually to 300 K in NVT ensemble using v-rescale algorithm for 10 ns. Then, NPT ensemble was used for 5 ns by restraining the complexes (unligated Mpro, Mpro–ECG, Mpro–EGCG, Mpro–GCG) while slowly allowing the solvent molecules to relax around it and finally another 10 ns NPT equilibration was carried out by gradually removing the restraints on the complexes. The pressure was maintained using Berendsen barostat. For every system, the average temperature and pressure values remained close to the desired values. The equilibrated systems were then subjected to unrestrained production MD simulations of 100 ns each, maintaining target pressure (1 bar) and temperature (300 K). Reproducibility of the results was checked by five different repeating simulations with different initial velocities and equilibration times. The overall trends in the RMSD, root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA) and hydrogen bond numbers remained the same. RMSD, RMSF, Rg, SASA and total number of H-bond count were calculated from the MD simulations (Purohit, 2014; Rajendran, 2016; Rajendran et al., 2018; Sharma et al., 2020; Singh et al., 2020).

MM–GBSA

To predict the theoretical free energies of binding of ligands to receptor, generally two most commonly used methods are (a) the molecular mechanics generalized Born surface area (MM-GBSA) and (b) molecular mechanics Poisson–Boltzmann surface area (MM-PBSA). Both methods are equally efficient to predict the correct binding affinities (Chen, 2016; Chen et al., 2015; Hou et al., 2011; Venugopal et al., 2020). Here, we used the MM-GBSA method to calculate the relative binding free energies of green tea polyphenols to Mpro protein. The free energy of binding can be calculated as ΔGbind = ΔH − TΔS.

ΔH = ΔEelec + ΔEvdW + ΔGpolar + ΔGnon-polar, where Eelec and EvdW are the electrostatic and van der Waal’s contributions, and Gploar and Gnon-polar are the polar and non-polar solvation terms, respectively. The non-polar energy is calculated from the SASA while the polar contribution of the free energy can be estimated using generalized Born model with an external dielectric constant of 80 and internal dielectric constant of 1. The entropic contribution is neglected here due to similar type of ligands bind to the receptor. Therefore, our calculated values are referred to as relative binding free energies. The theoretical binding free energy of the potent inhibitors of Mpro was identified using the prime module of Schrodinger suit (Schrödinger Release 2020-1: Prime, Schrödinger, LLC, New York, NY, 2020). MM-GBSA is a popular method to calculate binding energy, which uses energy properties of free ligand, free receptor and receptor–ligand complex for binding affinity calculation. Binding energies were estimated for the three polyphenols selected based on the binding affinity of AutoDock Vina docking and three complex structures were selected. Then with these docking structures MM-GBSA were calculated.

Pharmacokinetic properties analysis

Various pharmacokinetic properties of these eight polyphenols and N3 (chosen as model inhibitor) were predicted using pkCSM-pharmacokinetics and Swiss ADME online softwares (Daina et al., 2017; Pires et al., 2015). The drug-likeness properties such as absorption, distribution, metabolism and excretion parameters of these polyphenols along with the toxicity were mainly analyzed.

Result and discussion

The theoretical study of SARS CoV-2 found a new horizon when Yang and coworkers solved the crystal structure of Mpro with N3 complex (Jin et al., 2020). The crystal structure revealed that Mpro is a cysteine protease which is a homodimer consisting of two identical protomers. Each protomer has three distinct domains (domain I, II and III) along with a specific substrate-binding site as shown in Figure 1. This substrate-binding site mainly consists of a cysteine–histidine dyad (His41 and Cys145) which controls the catalytic activity of SARS CoV-2 Mpro. The substrate-binding pocket is located between domains I (magenta) and domain II (cyan) of each protomer (Figure 1). The compound N3 is reported to be a potential inhibitor of Mpro which mainly binds to the active/catalytic site of the protease Mpro. Docking of N3 compound with Mpro protease revealed that the structural orientation of N3 helps its binding to the catalytic cleft or substrate-binding region of Mpro (Figure 1 and Table 1).

Figure 1.

Surface representation of SARS CoV-2 Mpro with N3 inhibitor. The protomer of Mpro from SARS CoV-2 has been shown with its domains – Domain I (colored with magenta), Domain II (colored with cyan) and Domain III (colored with yellow) while blue color represents the linker. The inhibitor N3 (represented by the red stick) is attached to the substrate-binding pocket of Mpro.

Table 1.

Structure and binding energy of green tea polyphenols with Mpro along with N3 inhibitor as standard.

Complex	Structure	Binding energy (kcal/mol)
N3		−7.0
Epigallocatechin (EGC)		−7.0
Gallocatechin (GC)		−7.1
Catechin (C)		−7.1
Epicatechin (EC)		−7.2
Catechin gallate (CG)		−7.2
Epigallocatechin gallate (EGCG)		−7.6
Epicatechingallate (ECG)		−8.2
Gallocatechin-3-gallate (GCG)		−9.0

Molecular docking studies

Once the crystal structure of N3 inhibitor bound to Mpro was elucidated, many such theoretical studies commenced in order to find other suitable inhibitors for Mpro to control COVID-19 (Elmezayen et al., 2020; Huynh et al., 2020; Islam et al., 2020; Mahanta et al., 2020; Mittal et al., 2020). In the current study, we used this N3 complex as a standard inhibitor of Mpro and docked this complex with the protease. The binding energy of this particular binding was found to be −7.0 kcal/mol which is identical to the previously published estimation by other investigators (Huynh et al., 2020). Analysis of N3 docked complex suggested that this complex is stabilized by multiple hydrophobic and hydrogen bond interactions (particularly with His41 and Cys145) (Tables 2 and 3). Thus, we can also conclude like others that binding of N3 to both the catalytic residues (His41 and Cys145) of Mpro may be an important factor behind the inhibition of its protease activity. But the pharmacokinetics analysis revealed that N3 shows hepatotoxicity which makes it carcinogenic for humans (Table 4). Data presented in this table also reflected a negative tolerance dose of this compound for human. All these adverse pharmacokinetic behaviors make N3 an unsafe drug for COVID-19 treatment.

Table 2.

Hydrogen bond interactions of N3 and different green tea polyphenols with the SARS CoV-2 Mpro.

Complex	Number of H-bonds	Amino acids of Mpro i nvolved in H-bonding	Hydrogen bond distance (Å)
N3	8	His41	2.3
		Glu166	2.8
		Cys145	2.1
		Phe140	3.2
		Thr190	2.8
		His164	2.0
		Gly143	2.9
		Gln189	2.9
Epigallocatechin (EGC)	4	Ser144	2.5
		His163	2.9, 3.3
		Gln192	2.6
Gallocatechin (GC)	4	Phe140	2.4
		Glu166	2.3
		Arg 188	2.0
		Gln192	2.7
Catechin (C)	5	Leu141	2.
		Ser144	2.5
		His163	3.0, 3.3
		Gln192	2.6
Epicatechin (EC)	5	Ser 144	2.4
		His163	2.9, 3.3
		Thr190	2.2
		Gln192	2.6
Catechin gallate (CG)	5	Ser144	2.4
		His163	3.1
		Arg188	2.8
		Thr190	2.5, 1.9
Epigallocatechin gallate (EGCG)	9	Thr26	2.2, 1.9
		His41	2.8
		Cys145	2.6
		Ser144	2.3, 2.7
		Glu166	2.9
		Gln189	2.1
		Gly143	2.7
Epicatechingallate (ECG)	8	Ser144	2.4, 2.7
		Gly143	2.1
		Cys145	2.1
		Thr26	2.3, 1.9
		His41	2.8
		Glu166	2.9
Gallocatechin-3-gallate (GCG)	9	Phe140	2.1, 2.1
		His163	2.9
		Ser144	2.2, 2.5, 2.7
		Cys145	2.7
		Gly143	2.3, 2.7

Table 3.

Binding interactions of N3 and different polyphenols of green tea with the active site of SARS CoV-2 Mpro.

Complex	Non-covalent interactions (other than H-bonding)
	Polar	Hydrophobic
N3	Ser144, Asn142,	Phe140, Leu141, Met165, Glu166, Leu167, Pro168, His172, Ala191, Gln192
Epigallocatechin (EGC)	Ser46, Asn142	Leu141, Met165, Glu166, Ala191
Gallocatechin (GC)	Ser46, Asn142	Met49, Leu141, Met165, Gln189
Catechin (C)	Thr24, Thr25, Thr45, Ser46	Leu27, Met49
Epicatechin (EC)	Thr45, Ser46, Gln189	Met49, Leu141, Met165
Catechin gallate (CG)	Gln189	Phe140, Leu141, Met165, Glu166, His172
Epigallocatechin gallate (EGCG)	Thr24, Thr25, Ser46, Asn142, Gln192	Met49, Leu141, His163, Met165
Epicatechingallate (ECG)	Thr24, Thr25, Thr45, Ser46, Asn142, Gln189	Met49, Leu27, Leu141, Met165
Gallocatechin-3-gallate (GCG)	Asn142, Gln189, Thr190, Gln192	Met49, Leu141, Met165, Glu166, Arg188

Table 4.

Pharmacokinetic properties of N3 and green tea polyphenols.

Compound	MW	H-Ac	H-Do	Nrot	TPSA	LogP	IA	TC	LD50	HT	AT	MTD	NLV
N3	684.82	10	5	23	196.88	2.1352	62.162	0.843	4.355	Yes	No	−0.341	2
Epigallocatechin (EGC)	306.27	7	6	1	130.61	1.2517	54.128	0.328	2.492	No	No	0.506	1
Gallocatechin (GC)	306.27	7	6	1	130.61	1.2517	54.128	0.328	2.492	No	No	0.506	1
Catechin (C)	290.27	6	5	1	110.38	1.5461	68.829	0.183	2.428	No	No	0.438	0
Epicatechin (EC)	290.27	6	5	1	110.38	1.5461	68.829	0.183	2.428	No	No	0.438	0
Catechin gallate (CG)	442.37	10	7	4	177.14	2.5276	62.096	−0.169	2.558	No	No	0.449	1
Epigallocatechin gallate (EGCG)	458.37	11	8	4	197.37	2.2332	47.395	0.292	2.522	No	No	0.441	2
Epicatechin gallate (ECG)	442.37	10	7	4	177.14	2.5276	62.096	−0.169	2.558	No	No	0.449	1
Gallocatechin-3-gallate (GCG)	458.37	11	8	4	197.37	2.2332	47.395	0.292	2.522	No	No	0.441	2

MW, molecular weight (g/mol); H-Ac, no. of hydrogen bond acceptor; H-Do, no. of hydrogen bond donors; Nrot, no. of rotatable bonds; TPSA, topological polar surface area (Å2); LogP, predicted octanol/water partition coefficient; IA, intestinal absorption (% absorbed); TC, total clearance (log ml/min/kg); LD50, oral rat acute toxicity; HT, hepatotoxicity; AT, AMES toxicity; MTD, maximum tolerated dose for human (log mg/kg/day); NLV, no. of Lipinski rule violation.

In search of effective drugs for COVID-19 treatment, we hypothesized that green tea catechins or polyphenols (structure given in Table 1) may be suitable drug candidates by exerting inhibitory effects like N3. Their antiviral activity against diverse viruses including RNA viruses further strengthens our hypothesis (Calland et al., 2012; Carneiro et al., 2016; Ide et al., 2016; Ismail & Jusoh, 2017; Weber et al., 2015). But, before testing our hypothesis, first, we analyzed the pharmacokinetic and the drug-likeness properties of these eight green tea polyphenols. All the polyphenols exhibited negative AMES toxicity and hepatotoxicity which suggested that these compounds are not carcinogenic in nature (Table 4). Besides these, the oral rat acute toxicity (LD50) also supports this phenomenon (Table 4). All the polyphenols show good absorption in the human intestine as well as good excretion out (except catechin gallate and epicatechin gallate) from the body. Also, these polyphenolic compounds have a maximum tolerance dose (expressed in terms of log mg/kg/day) ranging from 0.438 to 0.506 which makes them suitable for human use. Altogether this valuable information indicated that green tea polyphenols/catechins can be safely used as drugs.

Once we understood that these polyphenols are safe, we studied their inhibition potency against SARS CoV-2 Mpro using molecular docking study. Amongst them, five polyphenols (EGC, GC, C, EC and CG) bound to Mpro having binding energy comparable to that of ‘Mpro–N3 binding’ (Table 1). However, these four polyphenols fitted within the substrate-binding pocket of Mpro having slightly different binding mode(s). For EGC–Mpro complex, two polar interactions (with Ser46 and Asn142 of Mpro), four hydrophobic interactions (with Leu141, Met165, Glu166 and Ala191 of Mpro) and four hydrogen bonds [with Ser144 (2.5 Å), His163 (2.9 Å & 3.3 Å) and Gln192 (2.6 Å) of Mpro] were observed (Supplemental Figure 1, Tables 2 and 3). When GC was docked into the active site of Mpro, we observed a similar number of polar, hydrophobic and hydrogen bond interactions, that is, two polar interactions (with Ser46, and Asn142 of Mpro), four hydrophobic interactions (with Met49, Leu141, Met165 and Gln189 of Mpro) and four hydrogen bonds [with Phe140 (2.4 Å), Glu166 (2.3 Å), Arg188 (2.0 Å) and Gln192 (2.7 Å) of Mpro] (Supplemental Figure 2, Tables 2 and 3). When the other three polyphenols (C, EC and CG) were docked individually to Mpro, these complexes were stabilized by more number of hydrogen bond interactions (5 nos.) and many non-covalent (polar and hydrophobic) interactions (Supplemental Figures 3–5, Tables 2 and 3). Leu141, Ser144, His163 and Gln192 are the four amino acid residues of Mpro involved in hydrogen bonding with C (Supplemental Figure 3, Table 2). However, EC formed hydrogen bonds with Ser 144, His163, Thr190 and Gln192 of Mpro and the amino acid residues of Mpro engaged in the hydrogen bonding with CG are Ser144, His163, Arg188 and Thr190, respectively (Supplemental Figures 4 and 5, Table 2). Critical analysis of docking results revealed that none of these five polyphenols have any interaction with the two key catalytic residues of Mpro (His41 and Cys145). Therefore, it can be concluded that these five polyphenols (EGC, GC, C, EC and CG) may not inhibit the catalytic/proteolytic activity of Mpro and thus may be ineffective towards the treatment of COVID-19 disease.

On the contrary, three other polyphenols (EGCG, ECG and GCG) interacted with one or both of the catalytic residues (His41 and Cys145) of Mpro by hydrogen bonding (Figures 2–4, Table 2). Both EGCG and ECG formed hydrogen bond with His41 and Cys145 amino acid residues of Mpro (Figures 2 and 3, Table 2). Besides these two key residues, several other amino acid residues from the active site of Mpro were involved in hydrogen bonding as well as different non-covalent interactions (Figures 2 and 3, Tables 2 and 3). Additionally, we noticed that ‘Mpro–GCG docked complex’ was stabilized by nine hydrogen bond interactions [with Phe140 (2.1 Å & 2.1 Å), His163 (2.9 Å), Ser144 (2.2 Å, 2.5 Å & 2.7 Å), Cys145 (2.7 Å) and Gly143 (2.3 Å & 2.7 Å) of Mpro] and multiple non-covalent interactions (polar, hydrophobic, etc.) (Figure 4, Tables 2 and 3). The binding energy of these three polyphenols are higher than that of ‘Mpro–N3 docked complex’ (−7.6 to −9.0 kcal/mol) (Table 1). Since, these three polyphenols (EGCG, ECG and GCG) interacted with the catalytic residue(s) of the Mpro protease and exhibited higher binding affinity than a well-known Mpro protease inhibitor (N3); we can say that they may possibly inhibit the catalytic activity of Mpro protease and maybe the good drug candidates for the treatment of COVID-19 disease.

Figure 2.

Molecular docking of EGCG with Mpro. Stereoview of the docked conformation of the Mpro–EGCG complex showing the possibility of hydrogen bonding interactions with the amino acid residues of Mpro (panel A). Surface representation showing the interaction of EGCG (green stick) at the substrate-binding region of Mpro (panel B). EGCG forms hydrogen bonding with many amino acid residues including His41 and Cys145 of Mpro.

Figure 3.

Molecular docking of ECG with Mpro. The docked conformation of the Mpro–ECG complex depicting the possible hydrogen bonding interactions with various amino acids of Mpro (panel A). Surface representation showing the binding of ECG (red) with Mpro (panel B). ECG forms hydrogen bonding with many amino acid residues including His41 and Cys141 of Mpro.

Figure 4.

Molecular docking of GCG with Mpro. Various hydrogen bonds with different amino acid residues of Mpro with GCG are shown in panel A as docked stereoview conformation. Binding of GCG (orange) at the active site of Mpro is illustrated in panel B as surface representation. GCG interacts with nine amino acid residues including Cys145 of Mpro via H-bonds.

In order to consolidate this conclusion, we selected Mpro–EGCG, Mpro–ECG and Mpro–GCG complex for MD study.

The MD simulations for the Mpro–EGCG, Mpro–ECG and Mpro–GCG complexes along with unligated Mpro/Mpro (unligated) and Mpro–N3 complex were performed for 100 ns. MD trajectories were initially analyzed with the aid of RMSD and RMSF so as to understand the stability and the fluctuations of these Mpro–polyphenol complex structures. The RMSD values of Mpro (unligated) and Mpro–ECG complex from 2 to 30 ns remained almost constant (∼1.3 Å) (Figure 5). After that RMSD values of these two systems increased gradually till 60 ns and thereafter it converged to a constant value (Figure 5). On the other hand, the RMSD values of Mpro–N3, Mpro–EGCG and Mpro–GCG complexes at 2 ns were 1.50, 1.55 and 1.45, respectively and these values remained almost the same till 30 ns. Afterward, RMSD values for these three complexes experienced a slight deviation (decrease for two Mpro–polyphenol complexes and increase for Mpro–N3 complex) up to 45 ns and beyond this time period, it converged to a constant value (Figure 5). The average RMSD values of Mpro (unligated), Mpro–N3, Mpro–ECG, Mpro–EGCG, Mpro–GCG complexes are 1.68, 1.56, 1.61, 1.53, 1.43 Å, respectively. Such low RMSD value (<2 Å) clearly suggested that all these three Mpro–polyphenol complexes are stable. From these data, it can be further inferred that the Mpro–GCG complex is more stable than the other two Mpro–polyphenol complexes. These findings additionally indicated that the stability of these three Mpro–polyphenol complexes are comparable or relatively more than that of Mpro–N3 complex.

Figure 5.

Determination of RMSD of unligated Mpro and Mpro–N3/ECG/EGCG/GCG complex. The MD simulations for each system were performed for 100 ns. These MD trajectories were analyzed with the aid of RMSD.

The RMSF of alpha-carbon atoms of all system were analyzed and represented in Figure 6. All the five systems [Mpro (unligated), Mpro–N3 and three Mpro–polyphenol complexes] exhibited a similar kind of fluctuation pattern across all three domains (Figure 6, Supplemental Figures 6–8). In all these systems, the fluctuations for the amino acid residues corresponding to domain I was lowest (Figure 6, Supplemental Figure 6). Whereas, the same for the amino acid residues corresponding to domain III was highest (Figure 6, Supplemental Figure 8). The average RMSF values of Mpro (unligated), Mpro–N3, Mpro–ECG, Mpro–EGCG, Mpro–GCG complexes are 1.28, 1.15, 1.22, 1.15, 1.17 Å, respectively. These values indicated that all the Mpro–polyphenol complexes as well as Mpro–N3 complex experienced relatively less conformational fluctuation than Mpro (unligated) system. The fluctuation of various specific amino acid residues (Thr25, Thr26, Leu27, His41 in domain I and Phe140, Leu141, Gly143, Ser144, Cys145, His163, His164, Glu166 of domain II) in these three Mpro–polyphenol systems were less than that in Mpro (unligated) system. Such less fluctuations of these residues in these three Mpro–polyphenol complexes further indicated that these residues within the active site of Mpro interact with ECG/EGCG/GCG. Moreover, these findings suggested that the two Mpro–polyphenol complexes (Mpro–EGCG and Mpro–GCG) and Mpro–N3 complex experienced comparable conformational fluctuation.

Figure 6.

Determination of RMSF of unligated Mpro and Mpro–N3/ECG/EGCG/GCG complex. The RMSF values for Mpro (unligated) and Mpro–N3/ECG/EGCG/GCG complex were estimated from the respective 100 ns MD trajectories. The values were then plotted separately for the domain I (amino acid residues 8–101) as shown in panel A, domain II (amino acid residues 102–184) as shown in panel B and domain III (amino acid residues 201–303) as shown in panel C.

MD trajectories corresponding to three Mpro–polyphenol and Mpro (unligated) systems were further analyzed with the aid of Rg, SASA and the existence of total number of intermolecular hydrogen bonds. Rg, which indicates the compactness of the system, was determined (Supplemental Figure 9(A)) and the calculated average Rg values are listed in Table 5. The obtained average Rg value for Mpro (unligated) was 2.5 Å. Interestingly, the Rg value for all these three Mpro–polyphenol complexes were slightly lower (2.20–2.24 Å) than that of Mpro (unligated) system. Therefore, it can be suggested that these Mpro–polyphenol complexes are relatively more rigid than the Mpro (unligated) system. Subsequently, the degree of expansion of protein volume in each system was assessed by estimating the SASA from individual MD trajectories (Supplemental Figure 9(B)) and the determined average SASA values are listed in Table 5. The SASA value corresponding to the ‘Mpro-ECG complex’ (14208.34 Å2) was slightly higher than that of Mpro (unligated) system (14124.45 Å2), suggesting a slight expansion of Mpro during the interaction with ECG. Whereas, the obtained SASA value for the ‘Mpro–EGCG complex’ and ‘Mpro–GCG complex’ 14113.53 and 13954.64 Å2, respectively, were a little lower compared to that of Mpro (unligated) system (14124.45 Å2) which indicated that Mpro suffers less expansion upon binding with EGCG and GCG.

Table 5.

Average values of the Rg, SASA and the total number of intermolecular hydrogen bond formed for the simulated systems.

System	Rg (Å)	SASA (Å2)	Total no. of intermolecular hydrogen bonds formed
Mpro (unligated)	22.52	14124.45	518
Mpro–ECG	22.22	14208.34	520
Mpro–EGCG	22.45	14113.53	528
Mpro–GCG	22.07	13954.64	540

The conformational stability of these three Mpro–polyphenol complexes was analyzed by estimating the total number of intermolecular hydrogen bonds formed during the entire simulation time span (Table 5). The average number of intermolecular hydrogen bonds in Mpro (unligated) and Mpro–ECG complex was almost equal suggesting similar conformational stability in these two systems (Table 5). However, more number of intermolecular hydrogen bonds was evidenced in other two Mpro–polyphenol complexes. Between these two complexes, the lowest number of intermolecular hydrogen bonds (528) was observed when EGCG was complexed with Mpro. Whereas, in the Mpro–GCG complex, the average number of hydrogen bonds was found to be 540 suggesting the highest conformational stability of this complex. These findings are in agreement with the binding energy data obtained from molecular docking studies.

MM-GBSA

Drugs that showed a high binding score in AutoDock Vina are further selected to calculate binding energy using the MM-GBSA method. The MM-GBSA free energy values of ECG–Mpro, EGCG–Mpro and GCG–Mpro are calculated from the MD trajectories. The MM-GBSA free energy values of GCG–Mpro, EGCG–Mpro and ECG–Mpro complexes are −53.5, −48.9 and −43.56 kcal/mol, respectively. From the MM-GBSA values (Table 6) it is clear that GCG–Mpro complex shows higher free binding energy than other complexes which is also supported by our AutoDock Vina binding energy values. Higher number of hydrogen bonds and hydrophobic interactions may play the key component in MM-GBSA value for GCG–Mpro complex system. Even lower SASA value in case of GCG–Mpro complex also indicated that the non-polar residues are buried in the solvent which may favor the stability of the complex through the synergistic effects of hydrogen bonding and hydrophobic interactions.

Table 6.

MM-GBSA values of different Mpro–polyphenol complexes.

System	Binding free energy (kcal/mol)
Mpro–GCG	−53.54
Mpro–EGCG	−48.92
Mpro–ECG	−43.56

Conclusion

In summary, our molecular docking study revealed that five polyphenols [EGC, GC, C, EC and CG] which do not interact with the His41 and Cys145 from catalytic dyad of Mpro, has binding energy more than −7.0 kcal/mol. Other three polyphenols [EGCG, ECG and GCG] has the interaction with one or both of these residues. The binding of these three polyphenols ranges between −7.6 and −9.0 kcal/mol with the lowest affinity for GCG and the highest affinity for EGCG. RMSD, RMSF, Rg and SASA investigations strongly support these findings. Even, binding free energy estimations using the MM-GBSA method also reveal that Mpro–GCG complex (−53.54 kcal/mol) is relatively more stable than Mpro–ECG (−48.92 kcal/mol) and Mpro–EGCG complex (−43.56 kcal/mol). Altogether, our findings reveal that green tea catechins/polyphenols (especially EGCG, ECG and GCG) can be potent anti-COVID-19 drug candidates. Additionally, this study opens up futuristic testing (in vitro and in vivo) possibilities of these three green tea polyphenols against COVID-19.