In silico analysis of the interactions of certain flavonoids with the receptor-binding domain of 2019 novel coronavirus and cellular proteases and their pharmacokinetic properties.

Coronavirus Disease 2019 (COVID-19) has infected more than thirty five million people worldwide and caused nearly 1 million deaths as of October 2020. The microorganism causing COVID-19 was named as Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2 or 2019-nCoV). The aim of this study was to investigate the interactions of twenty-three phytochemicals belonging to different flavonoid subgroups with the receptor binding domain (RBD) of the spike glycoprotein of 2019-nCoV, and cellular proteases [transmembrane serine protease 2 (TMPRSS2), cathepsin B and L (CatB/L)]. The compounds interacted more strongly with CatB and CatL than with the other proteins. Van der Waals and hydrogen bonds played an important role in the receptor-ligand interactions. As a result of RBCI (relative binding capacity index) analysis conducted to rank flavonoids in terms of their interactions with the target proteins, (-)-epicatechin gallate interacted strongly with all the proteins studied. The results obtained from molecular dynamics and molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) methods also supported this data. According to Lipinski's rule of five, (-)-epicatechin gallate showed drug-likeness properties. Although this molecule is not capable of crossing the blood-brain barrier (BBB), it was concluded that (-)-epicatechin gallate can be evaluated as a candidate molecule in drug development studies against 2019-nCoV since it was not the substrate of P-gp (P-glycoprotein), did not inhibit any of the cytochrome Ps, and did not show AMES toxicity or hepatotoxicity on eukaryotic cells.

Introduction

Coronavirus Disease 2019 (COVID-19), which first appeared in the Wuhan region of China and then spread all over the world rapidly, has infected more than thirty five million people worldwide and caused nearly 1 million deaths as of October 2020 ( Worldometers.info, 2020; F. Wu et al., 2020b; Zhu et al., 2020 ). The microorganism causing COVID-19 was named as new type Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2 or 2019-nCoV) ( F. Wu et al., 2020b; Zhou et al., 2020; Zhu et al., 2020 ). Some outbreaks caused by SARS-CoV and MERS-CoV in various parts of the world over the past decade were recorded. 2019-nCoV was the third important member of the coronavirus family that caused severe respiratory disease and human deaths (Zaki et al., 2012; Zhong et al., 2003).

The SARS-CoV and 2019-nCoV classified in the sarbecovirus subgenus of Coronaviridae were called human SARS-related coronaviruses (SARSr-CoV). Genome analysis has shown that the genetic material of this virus is a single-stranded RNA molecule consisting of approximately 30.000 nucleotides. As a result of genomic analysis, the genes encoding the four basic proteins of 2019-nCoV were clarified as follows: nucleocapsid protein, envelope protein, membrane protein, and spike glycoprotein (A. Wu et al., 2020a). Spike glycoprotein is actually a type I glycoprotein that extends out of the surface of the virus and is the first component to come into contact with the host cell. Since the spike glycoprotein is the main component the virus uses to bind to receptors on the host cell surface, it is of great interest in the development of the treatment strategies of COVID-19 (Luan et al., 2020).

Angiotensin converting enzyme 2, also known as ACE2, binds to the receptor binding domain (RBD) of SARS-CoV. Thus, this enzyme acts as the receptor that the virus uses to enter the cell (F. Li et al., 2005; W. Li et al., 2003). ACE2 shows a wide distribution in many tissues in the body, especially intestine, kidney, testis, liver, lungs, and heart. ACE2 is known to play an important role in regulating blood pressure. In addition, depending on the control of blood pressure, it also has important physiological tasks in the regulation of kidney and heart functions (Anguiano et al., 2017). Since the RBD of 2019-nCoV interacted with human ACE2, many studies have shown that human ACE2 plays a critical role in the entry of 2019-nCoV into cells (Letko et al., 2020; Zhou et al., 2020).

The receptor-ligand relationship is vital in host specificity and the life cycle of the virus. Some researchers claim that 2019-nCoV was transferred from bats to humans (Zhou et al., 2020). However, what we know about the intermediate host, which was thought to play a role in this transfer process, is still limited. Some studies suggested that pangolin may have played an important role in the evolution of 2019-nCoV (Lam et al., 2020; Wong et al., 2020). As can be understood from these speculative findings, the questions of which animal was involved in the evolution of 2019-nCoV and which other animal species were infected by the virus are still unanswered. The interaction between spike glycoprotein and ACE2 has been clarified as a result of studies on what amino acids bound to ACE2 in the receptor-binding domain of 2019-nCoV (Andersen et al., 2020). According to these studies, Leu455, Phe486, Glu493, Ser494, Asp501, and Tyr505 were the amino acids that interacted with mammalian ACE2. Considering the interaction between the 2019-nCoV spike glycoprotein and ACE2, it was understood that both SARS-CoV and 2019-nCoV have the capacity to infect many mammalian species such as Chinese hamsters, pangolin, dogs, cats, etc. (Luan et al., 2020).

Besides ACE2, the transmembrane protease, serine 2 (TMPRSS2) and cathepsins (CatB and L) were also thought to mediate the entry of both SARS-CoV and SARS-CoV-2 into the host cell. TMPRSS2 is an androgen-responsive serine protease that cleaves the spike glycoprotein from the S1/S2 region, facilitating viral entry and activation (Hoffmann et al., 2020). CatL is an important lysosomal endopeptidase. It plays a critical role in initiating protein breakdown. There is evidence that this cysteine protease mediated the entry of SARS-CoV and related viruses into the host cell (Hoffmann et al., 2020; Huang et al., 2006; Sudhan & Siemann, 2015).

Plants produce a wide variety of compounds to perform their normal physiological functions. These include important phytochemicals such as tannins, lignans, phytoalexins, polyphenols, and flavonoids. Polyphenols and flavonoids make up the largest group of secondary metabolites and are abundant in many sources such as spices, nuts, seeds, fruits, vegetables, stems, red wine, and tea that people consume frequently (Middleton et al., 2000). Plants synthesize these compounds in response to abiotic and biotic stress conditions such as UV rays, pathogens and insects (Carletti et al., 2014; Mandal et al., 2010). Studies have revealed that flavonoids have a wide range of biological/pharmacological properties such as antifungal, antibacterial, anticancer, anti-inflammatory, and antioxidant. In recent in vivo and in vitro studies, some flavonoids have been also found to exhibit considerable antiviral activity (Zakaryan et al., 2017).

The aim of this study was to determine the interaction of certain flavonoids (Figure 1) with the spike glycoprotein of 2019-nCoV, and host proteases (TMPRSS2, CatB and CatL). As a result of the receptor-ligand interactions, binding affinities and inhibition constant (Ki) values were determined. In order to rank flavonoids in terms of their interactions with all target proteins, RBCI (relative binding capacity index) analysis were conducted and ‘hit’ flavonoids were determined. Results of the docking analysis of the best molecules were also confirmed by subjecting them to the molecular dynamics simulations for 150 ns with all receptor proteins with subsequent estimation of the binding free energy using molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) methods. Drug-likeness properties according to the Lipinski's rule of five, ADMET (absorption, distribution, metabolism, excretion, and toxicity) profiles and intracellular target predictions of the flavonoids were also given.

Chemical structures of the flavonoids

Material and methods

Structural optimization of ligands

The protein data bank (pdb) files of all the ligands (quercetin, kaempferol, isorhamnetin, rutin, myricetin, resveratrol, naringenin, hesperidin, naringin, eriodictyol, genistein, daidzein, glycitein, apigenin, luteolin, tangeretin, epicatechin gallate, catechin, delphinidin, malvidin, peonidin, pelargonidin, petunidin) have been downloaded from PubChem (https://pubchem.ncbi. nlm.nih.gov) using Vega ZZ 3.2.0.9 download module. In the Vega ZZ platform, the atom types and electrical charges of the ligands were optimized with MMFF94 force field and Gasteiger-Marsili parameters (Pedretti et al., 2004).

Energy minimization of 2019-nCoV ACE2-RBD, TMPRSS2, CatB/L using nanoscale molecular dynamics (NAMD)

The structure of the spike glycoprotein was gained by removing the ACE2 subunit from the angiotensin-converting enzyme 2 - 2019nCoV RBD complex in the Vega ZZ environment. This model was downloaded from the following URL: https://swissmodel.expasy.org/interactive/HLkhkP/models/03 (PDBID: model_03. pdb) (Camacho et al., 2009; Remmert et al., 2012). Since the structure of the spike glycoprotein in model_03 shows a sequence identity of 100% to the 2019-nCoV ACE2 binding domain, this model was chosen as an appropriate 3 D structure in the molecular docking analyses. During the protein preparation step, the atom types and electrical charges of the spike glycoprotein were fixed using CHARMM22_PROT force field and Gasteiger-Marsili charges. Next, for the energy minimization of the spike glycoprotein with NAMD, each parameter was loaded from a template file. The number of time steps (number of minimization steps) were set to 10.000 and CHARMM22_PROT was set as the force field. When the energy minimization was completed, the 3 D structure corresponding to the last minimization step was saved as the lowest energy conformation. Also, to keep the spike glycoprotein structurally closer to the original crystallographic data, atom constraints were also applied to the protein backbone. In the energy minimization of TMPRSS2, CatB, and CatL, the same steps described above for the minimization of the 2019-nCoV spike glycoprotein were applied.

Homology modeling of TMPRSS2

The crystallographic data of TMPRSS2 enzyme structure has not been resolved until today, therefore we generated a homology model of this enzyme to utilize in further molecular docking and molecular dynamics analyses. The amino acid sequence of TMPRSS2 was downloaded from UniProtKB (https://www.uniprot.org/uniprot/O15393). Template search for TMPRSS2 catalytic domain was performed against the SWISS-MODEL template library with BLAST and HHBlits. BLAST was used to search the TMPRSS2 catalytic domain target sequence against the primary amino acid sequence in the SMTL. As a result of the BLAST search, a total of 788 templates were found. An initial HHblits profile has been built using the procedure as described in (Remmert et al., 2012). This procedure was followed by 1 iteration of HHblits against NR20. The obtained profile was then searched against all profiles of the SMTL and, finally, a total of 1167 templates were found.

ProMod3 was used to carry out model building for TMPRSS2 catalytic domain based on the target-template alignment. The rest of the procedure was carried out as previously described (Guex et al., 2009).

The model quality (global and per-residue) of TMPRSS2 obtained was evaluated with the QMEAN scoring function (Studer et al., 2020). A near-zero QMEAN score is a good value in terms of the quality of the fit between model structure and the experimental structure. According to the QMEAN score, however, scores of 4.0 and below indicate that the model is of poor quality. Therefore, among the top 5 TMPRSS2 models we obtained as a result of homology modeling, we determined the 5ce1.1.A (model 06) model as the target structure in the molecular docking analysis.

In addition, whether our model has an energetically favorable conformation, we generated a Ramachandran plot (Figure S2) using the PROCHECK web server (Laskowski et al., 1993). Also, ERRAT web-based tool (Figure S1) was also deployed to calculate the overall quality factor (OQF) for non-bonded atomic interactions (Colovos & Yeates, 1993).

Figure 2.

RBCI values of the flavonoids.

Molecular docking analyses

Firstly, we thought that it could be helpful to explain more clearly the expressions of RBM (receptor binding motif) and RBD (receptor binding domain) in the literature: Spike's RBD is the smallest part of this protein which is supposed to be able to capture the interaction with the receptor, hACE2. The RBD itself is on its hand composed by a core (residues 333–438 and residues 507–527) and a Receptor Binding Motif (RBM, residues 438–506). This RBM is the part of the RBD which is directly interacting with the hACE2 ( de Andrade et al., 2020; Lan et al., 2020; Omotuyi et al., 2020; Y. Wang et al., 2020; Woo et al., 2020 ). In this sense, throughout this study, we envisioned that RBM is a kind of "active site" of the RBD.

Molecular docking analyses between the target structures and the ligands were performed using AutoDock 4.2.6 and the corresponding docking scores (binding free energies) of the ligands with 2019-nCoV RBM, TMPRSS2, CatB (PDB ID: 1GMY), and CatL (PDB ID: 2YJ9) were calculated. AutoDockTools-1.5.6 was used to prepare the target and ligand molecules and also the parameters prior to initiating the docking analysis using AutoDock 4.2.6 (Sanner, 1999). In this study, the grid box coordinates used in molecular docking analyzes were adjusted to ensure that all the tested phytochemicals interact with amino acids in the active sites of the enzymes in question (Andersen et al., 2020; Greenspan et al., 2001; Hardegger et al., 2011; Wilson et al., 2005).

Prior to molecular docking analyzes, polar hydrogen atoms in the receptor and the ligand molecules were retained, while nonpolar hydrogens were merged and then, the Gasteiger charges of the ligands were calculated with AutoDockTools (Morris et al., 2009; Nasab et al., 2017; Sanner, 1999). Also, the Kollmann charges were added for the receptor. During the docking experiments, all the rotatable bonds of the ligands were allowed to rotate and then the optimized protein (rigid) and ligand (flexible) structures were saved in PDBQT format. Grid box coordinates were adjusted as: a) 80 × 90 × 40 Å points for the spike glycoprotein; b) 60 × 110 × 86 Å points for TMPRSS2; c) 86 × 84 × 44 Å points for CatB; and d) 54 × 52 × 60 Å points for CatL. These grid box sizes were previously determined to cover the active amino acid residues of the enzymes in question.

In all docking analyses, 50 genetic algorithm (GA) runs using an initial population of 150 individuals, maximum number of 2.500.000 energy evaluations, and a maximum number of 27.000 generations were selected. The values of 0.02 and 0.8 were chosen as the default parameters for mutation and crossover rates, respectively. After 50 independent docking runs, all the possible binding modes (conformations) of the ligands were clustered by the program and were ranked based on the most negative binding free energy (kcal/mol) of the ligand conformation. The best docking poses obtained using the AutoDock 4.2.6 between the ligand and receptor structures were analyzed with the BIOVIA Discovery Studio Visualizer 2016.

Furthermore, according to the results obtained from molecular docking experiments, heatmaps of ligands (Supplementary Figures S3, S4, S5 and S6) that interact with specific amino acid residues of each receptor were created. Heatmaps are built on the logic of how many times each ligand interacts with each residue (regardless of bond type). Thus, heatmaps provide an overview of the frequency of interactions of 23 different ligands with certain amino acids of these receptors.

Success criteria set in docking analysis

In the current study, the lowest of all clusters in terms of the binding free energy was considered as the energetically most stable configuration and taken as reference (Morris & Lim-Wilby, 2008). The calculated inhibition constant (Ki) obtained with AutoDock 4.2.6 for each docked phytochemical were also represented.

Calculation of relative binding capacity index (RBCI) values

In this study, we have developed a new analysis method called RBCI and applied it to statistically rank the activity potentials of phytochemicals by using binding free energy and inhibition constant values obtained from the docking analyses (Istifli et al., 2020). Using the RBCI parameter, it is possible to compare statistically related data with different scientific meanings. If the ranking of the interactions of the ligands with proteins is done only in light of one parameter (e.g. binding free energy or inhibition constant only), the molecules can only be ordered in terms of their potential in this parameter. However, ranking based on only one of these parameters cannot represent the full activity potential of these molecules. The most common method used to calculate the interaction between the receptor and the ligand in multiple measurements is the 'central tendency', in which the components are ranked based on the average value of each component. However, since the units and scales of the data obtained from each parameter are different, it is not possible to obtain a universal value for all components.

If the values in each data set (binding energy and inhibition constant) are converted to standard scores, it is possible to compare them with each other. In order to calculate the arithmetic mean values, first of all, binding energy and inhibition constant data of each phytochemical were used regardless of their units and raw values were obtained. These raw values calculated for each component were subtracted from the arithmetic mean and divided by standard deviation, and standard scores were obtained (see equation given below) (Sharma, 1996). RBCI values of each phytochemical were calculated by averaging these standard scores obtained separately for each protein target.where ‘x’ is the raw data, ‘μ’ is the mean, and ‘σ’ is the standard deviation.

Although RBCI is a relative index and does not represent the specific binding capacities of the components, it makes it possible to rank components reasonably based on their binding energy and inhibition constant values. Therefore, it can be used as an integrated approach to evaluate the molecular interaction of the components, considering all parameters.

Molecular dynamics of top-ranked receptor-ligand complexes

The configurations with the highest affinity (most negative binding free energy) obtained from the docking of receptor-ligand complexes were stored and prepared for the molecular dynamics simulations as follows: The receptor coordinates files were converted from pdbqt format to pdb format using Babel (O’Boyle et al., 2011) and thereafter to use in molecular dynamics simulations using standard GROMACS (Abraham et al., 2015) tools, choosing the AMBER03 force field (Duan et al., 2003). The ligand coordinates files were read with AutoDockTools (Morris et al., 2009) to add the hydrogen atoms, and the obtained pdb files were then submitted to the Acpype online server (Sousa da Silva & Vranken, 2012) to generate GROMACS topology files, employing the parameters of the General Amber Force Field, GAFF (J. Wang et al., 2004) and AM1-BCC partial charges (Jakalian et al., 2002). The topologies and coordinates of receptors and ligands were then combined to construct the topologies and coordinates of the complexes. The complexes were then solvated using TIP3P water molecules (Jorgensen et al., 1983) as well as sodium and chloride ions corresponding to physiological concentration in cubic boxes under periodic boundary conditions. The van der Waals interactions were calculated until a cutoff of 1.2 nm and the electrostatic interactions were calculated using the PME method (Darden et al., 1993; Essmann et al., 1995). The systems were initially energetically minimized using the steepest descent algorithm, submitted to a 500 ps simulation employing position restraints for both receptor and ligand and to a sequence of three unrestricted molecular dynamics simulations with 5 ns each in temperatures of 200 K, 240 K and 280 K for the thermalization (heating) of the system. After these steps, 150 ns long production simulation runs were carried out. The simulation sampled the NPT ensemble, employing a Nosé-Hoover thermostat (Hoover, 1985; Nosé, 1984) and a Parrinello-Rahman barostat (Parrinello & Rahman, 1981). The trajectories were analyzed to quantify the structural and thermodynamic stability of the complexes and to identify the pattern of intermolecular interactions.

Detailed qualitative visual analyses of the time evolution of the interactions in the trajectories were also carried out. The strength of the interactions between ligands and receptors was estimated with Molecular mechanics Poisson-Boltzmann Surface Area (MM/PBSA) binding free energy calculations (Baker et al., 2001; Homeyer & Gohlke, 2012). Sets of 150 configurations for each system were obtained from 1 ns spaced snapshots obtained from the molecular dynamics trajectories. The calculations were carried out with g_mmpbsa (Kumari et al., 2014) using a gridspace of 0.5 A, salt concentration of 0.150 M, solute dielectric constant of 2 and using the solvent acessible surface area (SASA) as estimate of the nonpolar solvation energy.

Drug-likeness, ADMET profile and target prediction

The drug-likeness, ADMET and target profiles of potential hit compounds are very important in terms of reducing side effects in the pharmaceutical industry. In the current study, web-based SwissADME, SwissTargetPrediction and pkCSM tools were used to determine such effects of flavonoids analyzed (Daina et al., 2017, 2019; Pires et al., 2015).

Results and discussion

In this study, the molecular interactions of twenty-three phytochemicals belonging to four different flavonoid subgroups (flavonols, flavanones, isoflavones, and anthocyanidins) (Table 1, Figure 1) with the RBD of the spike glycoprotein of 2019-nCoV, TMPRSS2, CatB and CatL were analyzed.

Table 1.

PubChem CID, molecular weight and molecular formula of the compounds.

No	Compound	PubChem CID	Molecular weight (g/mol)	Molecular formula
Flavonols
1	Quercetin	5284452	302.23	C15H10O7
2	Kaempferol	5280863	286.24	C15H10O6
3	Isorhamnetin	44259381	478.40	C22H22O12
4	Rutin	5280805	136.23	C27H30O16
5	Myricetin	5281672	318.23	C15H10O8
6	Resveratrol	445154	228.24	C14H12O3
Flavanones
7	Naringenin	932	272.25	C15H12O5
8	Hesperidin	10621	610.60	C28H34O15
9	Naringin	442428	580.50	C27H32O14
10	Eriodictyol	102370911	450.40	C21H22O11
Isoflavones
11	Genistein	5280961	270.24	C15H10O5
12	Daidzein	5281708	254.24	C15H10O4
13	Glycitein	5317750	284.26	C16H12O15
14	Apigenin	5280443	270.24	C15H10O5
15	Luteolin	5280445	286.24	C15H10O6
16	Tangeretin	68077	372.40	C20H20O7
17	(-)-Epicatechin gallate	107905	442.40	C22H18O10
18	(+)-Catechin	9064	290.27	C15H14O6
Anthocyanidins
19	Delphinidin	68245	338.69	C15H11ClO7
20	Malvidin	159287	331.30	C17H15O7+
21	Peonidin	441773	301.27	C16H13O6+
22	Pelargonidin	440832	271.24	C15H11O5+
23	Petunidin	73386	352.72	C16H13ClO7

Molecular docking analysis

In receptor-ligand interactions, binding energy can vary depending on the three-dimensional structure of the receptor, the type of amino acids present in the active center and their locations. Therefore, different binding energy profiles were detected in each of the receptor-ligand pairs analyzed in the current study (Tables S1, S2, S3, and S4). According to these results, (-)-Epicatechin gallate, naringin, peonidin, pelargonidin showed high affinity to target proteins. On the other hand, rutin was one of the compounds with the lowest affinity for almost all protein targets. Ligands showed higher affinity to CatB and CatL than other proteins.

In order to determine the effective concentrations of ligands on protein targets, inhibition constants were calculated by AutoDock 4.2.6. The predicted inhibition constants of the flavonoids on the spike glycoprotein, TMPRSS2, CatB, and CatL were found to be in the ranges of 0.02-0.58, 0.04-339.73, 0.001-11.82 and 0.001-11.04 mM, respectively. (-)-Epicatechin gallate was found to be highly effective against all receptors.

Van der Waals interactions played an important role in receptor-ligand interplay. Classical and non-classical H bonds, hydrophobic and electrostatic interactions were also frequently found.

In the interaction between flavonoids and the RBD of the spike glycoprotein, the interaction with the amino acid residues of the receptor binding motif (RBM) inside the RBD such as Phe486 and Leu455 remained weak, while the ligands interacted with other residues of the RBM such as Tyr505, Asn501, Ser494, and Gln493 more effectively. Flavonoids also interacted intensely with Arg403, which is not one of the amino acid residues of the RBM of the spike glycoprotein (Figure S3).

Figure 3.

Top ranked conformations of (-)-epicatechin gallate. (A- RBM of the spike glycoprotein of 2019-nCoV, B- TMPRSS2, C- CatB, D- CatL).

Flavonoids interacted with His296 and Ser441, active amino acid residues of TMPRSS2, but did not establish any interactions with Asp 345. Other amino acids that interacted with flavonoids were Val280, Lys300, Tyr337, Lys340, Thr341, and Lys342 (Figure S4).

Figure 4.

Top ranked conformations of naringin. (A- RBM of the spike glycoprotein of 2019-nCoV, B- TMPRSS2, C- CatB, D- CatL)

Cys29, Gly27, His111, His199, and Trp221, active amino acid residues, played an important role in the interaction of flavonoids with CatB (Figure S5).

Figure 5.

Top ranked conformations of peonidin. (A- RBM of the spike glycoprotein of 2019-nCoV, B- TMPRSS2, C- CatB, D- CatL).

Cys25, Trp26, Gly68, Leu69, Met70, Ala135, Met161, and Asp162, active amino acid residues, played an important role in flavonoid/CatL interactions. Asp114, His163, Gly164, Ala214, Ala215, and Ser216, which are not active amino acid residues of CatL, also contributed significantly to these interactions (Figure S6).

Figure 6.

Top ranked conformations of pelargonidin. (A- RBM of the spike glycoprotein of 2019-nCoV, B- TMPRSS2, C- CatB, D- CatL).

As a result of the RBCI analysis, (-)-epicatechin gallate was the most effective molecule against all protein targets (Figure 2). It was followed by naringin, peonidine and pelargonidine, respectively. The top ranked conformations of these flavonoids (‘hit’ compounds) were given as Figures 3–6, respectively.

Drug-likeness properties and ADMET profiles of flavonoids

The drug-likeness properties of flavonoids were given in Table 2. Although rutin, hesperidin and naringin have three violations (MW > 500, N or O > 10, NH or OH > 5) and myricetin, (-)-epicatechin gallate and delphinidin have one violation (NH or OH > 5), all other flavonoids were found to have drug-likeness properties without any violation.

Table 2.

Drug-likeness properties of docked flavonoids.

No	Compound	Number of rotatable bonds	TPSA1	Consensus Log P	Log S (ESOL2)	Drug-likeness (Lipinski’s rule of five)
1	Quercetin	1	131.36 Å2	1.23	−3.16	Yes (0 Violation)
2	Kaempferol	1	111.13 Å2	1.58	−3.31	Yes (0 violation)
3	Isorhamnetin	2	120.36 Å2	1.65	−3.36	Yes (0 violation)
4	Rutin	6	269.43 Å2	−1.12	−3.30	No (3 violations; MW > 500, N or O > 10, NHorOH > 5)
5	Myricetin	6	151.59 Å2	0.79	−3.01	Yes (1 Violation; NH or OH > 5)
6	Resveratrol	2	60.69 Å2	2.48	−3.62	Yes (0 Violation)
7	Naringenin	1	86.99 Å2	1.84	−3.49	Yes (0 Violation)
8	Hesperidin	7	234.29 Å2	−0.72	−3.28	No (3 violations; MW > 500, N or O > 10, NH or OH > 5)
9	Naringin	6	225.06 Å2	−0.79	−2.98	No (3 violations; MW > 500, N or O > 10, NH or OH > 5)
10	Eriodictyol	6	107.22 Å2	1.45	−3.72	Yes (0 Violation)
11	Genistein	1	90.90 Å2	2.04	−3.72	Yes (0 Violation)
12	Daidzein	1	70.67 Å2	2.24	−3.53	Yes (0 Violation)
13	Glycitein	2	79.90 Å2	2.30	−3.57	Yes (0 Violation)
14	Apigenin	1	90.90 Å2	2.11	−3.71	Yes (0 Violation)
15	Luteolin	1	111.13 Å2	1.73	−3.71	Yes (0 Violation)
16	Tangeretin	6	76.36 Å2	3.02	−4.11	Yes (0 Violation)
17	(-)-Epicatechingallate	4	177.14 Å2	1.23	−3.70	Yes (1 violation; NH or OH > 5)
18	(+)-Catechin	1	110.38 Å2	0.85	−2.22	Yes (0 Violation)
19	Delphinidin	1	134.52 Å2	−0.79	−3.89	Yes (1 violation; NH or OH > 5)
20	Malvidin	3	112.52 Å2	0.71	−2.86	Yes (0 Violation)
21	Peonidin	2	103.29 Å2	0.76	−2.81	Yes (0 Violation)
22	Pelargonidin	1	94.06 Å2	0.73	−2.76	Yes (0 Violation)
23	Petunidin	2	123.52 Å2	−0.55	−3.36	Yes (0 Violation)

TPSA: Topological polar surface area.

ESOL: Estimated aqueous solubility.

Data source: https://www.swissadme.ch.

ADMET profiles of flavonoids were also given in Table 3. According to the data given in the table, it was understood that only quercetin, resveratrol, daidzein, and tangeretin could cross the blood-brain barrier. Rutin, hesperidin, naringin, (-)-epicatechin gallate, (+)-catechin, delphinidin, and petunidine did not have inhibitory effect on any of cytochrome P. None of the flavonoids showed hepatotoxicity. Apart from resveratrol, all other flavonoids were also found to show no AMES toxicity. LD50 values of the compounds in rats were in the range of 1.791-2.558 mol/kg.

Table 3.

ADMET profiles of flavonoids.

No	Compound	BBB1 permeation*	P-gp substrate2*	CYP3 inhibitionϯ	AMES Toxicityϯ	Hepatotoxicityϯ	LD50 in rat (mol/kg)ϯ
1	Quercetin	Yes	No	Yes (CYP1A2, CYP2D6, CYP3A4)	No	No	2.471
2	Kaempferol	No	No	Yes (CYP1A2, CYP2D6, CYP3A4)	No	No	2.449
3	Isorhamnetin	No	No	Yes (CYP1A2, CYP2D6, CYP3A4)	No	No	2.407
4	Rutin	No	Yes	No	No	No	2.491
5	Myricetin	No	No	Yes (CYP1A2, CYP3A4)	No	No	2.497
6	Resveratrol	Yes	No	Yes (CYP1A2, CYP2C9, CYP3A4)	Yes	No	2.529
7	Naringenin	No	Yes	Yes (CYP1A2, CYP3A4)	No	No	1.791
8	Hesperidin	No	Yes	No	No	No	2.506
9	Naringin	No	Yes	No	No	No	2.495
10	Eriodictyol	No	Yes	Yes (CYP3A4)	No	No	2.030
11	Genistein	No	No	Yes (CYP1A2, CYP2D6, CYP3A4)	No	No	2.268
12	Daidzein	Yes	No	Yes (CYP1A2, CYP2D6, CYP3A4)	No	No	2.164
13	Glycitein	No	No	Yes (CYP1A2, CYP2D6, CYP3A4)	No	No	2.170
14	Apigenin	No	No	Yes (CYP1A2, CYP2D6, CYP3A4)	No	No	2.450
15	Luteolin	No	No	Yes (CYP1A2, CYP2D6, CYP3A4)	No	No	2.455
16	Tangeretin	Yes	No	Yes (CYP2C9, CYP3A4)	No	No	2.368
17	(-)-Epicatechin gallate	No	No	No	No	No	2.558
18	(+)-Catechin	No	Yes	No	No	No	2.428
19	Delphinidin	No	Yes	No	No	No	2.548
20	Malvidin	No	Yes	Yes (CYP1A2)	No	No	2.346
21	Peonidin	No	Yes	Yes (CYP1A2)	No	No	2.408
22	Pelargonidin	No	Yes	Yes (CYP1A2)	No	No	2.432
23	Petunidin	No	Yes	No	No	No	2.453

BBB: Blood Brain Barrier.

P-gp: P-glycoprotein substrate.

CYP: Cytochrome P.

https://www.swissadme.ch.

http://biosig.unimelb.edu.au/pkcsm/prediction.

Target predictions of ‘hit’ flavonoids

The possible intracellular targets of (-)-epicatechin gallate, naringin, peonidine, and pelargonidine, were given in Figure 7(A–D), respectively. (-)-Epicatechin gallate appeared to act on proteases and kinases. While naringin was effective on proteases, oxidoreductases and electrochemical transporters, peonidine had an effect on lyase. Pelargonidine also appeared to be effective on many intracellular enzymes, especially oxidoreductases.

Figure 7.

Target predictions of A- (-)-epicatechin gallate, B- Naringin, C- Peonidin, and D- Pelargonidin.

Results of the molecular dynamics analysis of ‘hit’ flavonoids

The complex (-)-epicatechin gallate (ECG)-spike protein remained stable and the ECG molecule remained in the proximity of the receptor during all the simulation, with small rearrangements (Figure 8(A)). The ECG molecule started interacting with the RBM of the spike protein (mainly with residues Gln493, Ser494, Gly496, and Tyr505), but also with Glu406, Lys417 and Tyr453. During the simulation, the ECG molecule remained interacting in the RBM, which suffered only small rearrangements. At the end of the simulation the ECG molecule interacted with essentially the same residues as in the beginning, losing only the interactions with Lys417 and Gly496, but interacting with Ile418. The number of hydrogen bonds remained roughly conserved through the simulation (Figure 8(A)).

Figure 8.

Time evolution of the structure and interactions in the complexes of (-)-epicatechin gallate (ECG) with A- RBD of the spike glycoprotein of 2019-nCoV, B- TMPRSS2, C- CatB, and D- CatL in molecular dynamics simulations. In the left the root-mean-square deviation (RMSD) of atomic positions of the receptor (black) and of the complex ligand + receptor (red) are shown, whereas in the right the number of hydrogen bonds between ligand and receptor is shown.

The complex ECG-TMPRSS2 was stable during the simulation, with the ligand keeping essentially the same place interacting with the receptor (Figure 8(B)). The ECG molecule started interacting mainly with the residues His296, Tyr337, Asp338, Lys340, Thr341, Leu419 and Trp461 of the protein TMPRSS2. The molecule changed only slightly the position, although some structural rearrangements took place. At the end of the simulation the molecule interacted with the same residues as in the beginning (but losing the interaction with Asp338). There was some fluctuation on the number of hydrogen bonds along the simulation, nevertheless the interactions remained strong (Figure 8(B)).

The ECG molecule interacted strongly with the CatB, maintaining the distance and the majority of the interactions along the simulation (Figure 8(C)), as well as the high number of H-bonds. The ECG molecule remained almost in the same place: it started interacting with the residues Ser25, Cys26, Gly27, Asn72, Gly73, His110, Gly121, Glu122, Ser178, Met196, Gly197 and Gly198 of the receptor. At the end of the simulation, the interactions with Ser25 and Cys26 were lost, but new interactions with Trp30, His199 and Ala200 appeared, with the ECG molecule migrating slightly from the surface to a position more buried inside the protein. However, there were almost no structural changes on the protein (Figure 8(C)).

ECG and CatL, the interactions were moderately strong and some structural changes occurred (Figure 8(D)). There was some loss of hydrogen bonds along the simulation, with a recovery at the end. In the beginning of the simulation the ECG molecule was interacting mainly with the residues Ser4, Val5, Asp6, Trp7, Arg8, Glu9 and Lys10. The ECG molecule migrated at about 135 ns, losing contacts with the residues 4 until 9 and starting to interact with Gly11, Tyr12, Val13, Thr14 and Pro15 (Figure 8(D)).

The binding free energies calculated using the MM/PBSA method (Table 4) agree with the findings of the docking and molecular dynamics: there are favorable interactions between ECG and all receptors, as can be seen in the negative binding free energies. The strongest interaction was found for ECG with CatB and the weakest for ECG with CatL, in agreement with the time evolution of the number of hydrogen bonds and with the qualitative analyses of the simulations.

Table 4.

Binding free energies and their components, calculated using MM/PBSA, for the interaction between ECG and the receptors A- Spike, B- TMPRSS2, C- CatB and D- CatL.

Ligand	Protein	Van der Waals Energy	Electrostatic Energy	Polar Solvation Energy	SASA Energy	Binding Energy
ECG	Spike	−116.31 ± 2.30 kJ/mol	−81.18 ± 2.82 kJ/mol	160.83 ± 2.53 kJ/mol	14.71 ± 0.15 kJ/mol	-51.43 ± 1.38 kJ/mol
ECG	TMPRSS2	−121.51 ± 1.72 kJ/mol	−93.07 ± 2.89 kJ/mol	173.74 ± 3.06 kJ/mol	−14.87 ± 0.13 kJ/mol	-55.67 ± 1.60 kJ/mol
ECG	CatB	−173.72 ± 2.00 kJ/mol	−119.03 ± 1.71 kJ/mol	202.52 ± 2.32 kJ/mol	−19.02 ± 0.13 kJ/mol	-109.15 ± 1.48 kJ/mol
ECG	CatL	−88.08 ± 1.50 kJ/mol	−75.10 ± 3.21 kJ/mol	155.12 ± 4.02 kJ/mol	−12.13 ± 0.12 kJ/mol	-20.11 ± 1.99 kJ/mol

Results of the molecular dynamics and binding free energy analysis of other ‘hit’ flavonoids (naringin, peonidin and pelargonidin) were also given in Figures S7, S8 and S9 and Table S6, respectively in the Supplementary file.

It is impossible to discuss the literature data on all flavonoids studied in the current study. Instead, it would be more logical to discuss the literature data of ‘hit’ flavonoids on target proteins. Literature data obtained from the docking analysis of the phytochemicals analyzed in the current study was also presented in Table S5 in detail (Supplementary file).

Some recent studies have shown that catechins can be used as alternative agents in COVID-19 infection. According to Maiti and Banerjee (2020), catechins, such as epigallocatechin gallate or theaflavin gallate, interact more strongly with the amino acids of the spike protein than hydroxychloroquine. There are also some data in the literature that some catechins, such as theaflavin digallate, epigallocatechin gallate and gallocatechin gallate, inhibit the main protease of 2019-nCoV together with the spike glycoprotein (Manish, 2020; Peterson, 2020; Sayed et al., 2020; Tallei et al., 2020). In addition to 2019-nCoV, molecular docking studies performed on (-)-epicatechin 3-O-(3′-O-methyl) gallate, a green tea polyphenol, showed that the molecule in question can inhibit TMPRSS, a cellular protease (Rahman et al., 2020). Catechins such as (-)-catechin-3-gallate, (+)-epicatechin-3-gallate and (-)-epigallocatechin-gallate also have an inhibitory effect on cathepsins, a group of lysosomal enzymes (Majumdar et al., 2011; Zhang et al., 2012). In a study carried out by Devika and Prince (2008), a significant decrease in cathepsin B and D activities were detected in Wistar rats given epigallocatechin gallate at a dose of 30 mg/kg for 21 days compared to the control group. However, there are no reports in the literature directly about the effect of epicatechin gallate itself on CatB and especially CatL. Therefore, data presented in the current study on the interaction of (-)-epicatechin gallate with CatB and CatL could be assumed as the first record for the literature.

In the literature, there were no reports of the interaction of peonidine and pelargonidine with target proteins. On the other hand, only one study was found regarding the interaction of naringin with the RBD of the spike protein of 2019-nCoV (Bhowmik et al., 2020). According to the results of this study conducted using PyRx and iGEMDOCK, it was determined that the binding affinity of naringin to spike protein was −8.3 kcal/mol.

Conclusion

In this computational study, the interactions of 23 different phytochemicals belonging to different flavonoid subgroups with the RBD of the spike glycoprotein of 2019-nCoV, TMPRSS2, CatB and CatL were analyzed. Anthocyanidins, isoflavones and flavanones were found to have more powerful interactions with the target proteins. Considering the interactions of phytochemicals with all target proteins, it was concluded that (-)-epicatechin gallate was more effective than other compounds. This molecule showed drug-likeness properties according to Lipinski's rule of five. Although (-)-epicatechin gallate cannot cross the blood-brain barrier, it was concluded that it can be considered as a candidate molecule in the drug development processes against 2019-nCoV since it did not have any toxic effect on eukaryotic cells.

Click here for additional data file.