Theobroma cacao L. compounds: Theoretical study and molecular modeling as inhibitors of main SARS-CoV-2 protease.

Cocoa beans contain antioxidant molecules with the potential to inhibit type 2 coronavirus (SARS-CoV-2), which causes a severe acute respiratory syndrome (COVID-19). In particular, protease. Therefore, using in silico tests, 30 molecules obtained from cocoa were evaluated. Using molecular docking and quantum mechanics calculations, the chemical properties and binding efficiency of each ligand was evaluated, which allowed the selection of 5 compounds of this series. The ability of amentoflavone, isorhoifolin, nicotiflorin, naringin and rutin to bind to the main viral protease was studied by means of free energy calculations and structural analysis performed from molecular dynamics simulations of the enzyme/inhibitor complex. Isorhoifolin and rutin stand out, presenting a more negative binding ΔG than the reference inhibitor N-[(5-methylisoxazol-3-yl)carbonyl]alanyl-l-valyl-N~1~-((1R,2Z)-4-(benzyloxy)-4-oxo-1-{[(3R)-2-oxopyrrolidin-3-yl]methyl}but-2-enyl)-L-leucinamide (N3). These results are consistent with high affinities of these molecules for the major SARS-CoV-2. The results presented in this paper are a solid starting point for future in vitro and in vivo experiments aiming to validate these molecules and /or test similar substances as inhibitors of SARS-CoV-2 protease.

Introduction

The cocoa tree is known since ancient times in Mesoamerica, as drinks obtained from its seeds can be traced through archeological and chemical studies that show the presence, in pottery containers dating from approximately 1000 B.C., of compounds originating from the secondary metabolism of this species [1], [2]. Shortly after the arrival of the Spanish in the new world, this product became known in Europe, and in 1737 Carl Nilsson Linnaeus designated it taxonomically as Theobroma cacao L. [3]. Currently the genus Theobroma comprises 22 species and is classified within the family Malvaceae, and 9 of them are native to the Amazon forest [4], [5], [6]. Cocoa beans are a rich source of phenolic compounds (10–12%, dry weight), specifically catechins (flavan-3-ols) and procyanidins, which are a large and diverse group of metabolites acting as antioxidants [7], [8]. According to the literature reviewed, 2 main genetic groups are considered, known as Criollo and Forastero, based on their geographic origins and morphological characteristics. In addition, there is a third group known as Trinitario, which is considered a hybrid product of the previous varieties [9]. The phenol content depends mainly on the genetics of the trees and can also vary according to the area of cultivation, fruit maturity, climatic conditions, harvest time and storage time of the cobs until they are opened [10], [11], [12]. The 3 main groups of polyphenols are present in cocoa beans – catechins, anthocyanins and proanthocyanidins – representing 35%, 4% and 58%, respectively [13]. It has been shown that during the different stages of processing and manufacturing (fermentation, drying, roasting and refining), the chemical composition of cocoa beans changes drastically, affecting the content and structure of polyphenols, and resulting in a reduction of their antioxidant activity [2], [14], [15]. The high content of antioxidant compounds makes cocoa a very interesting product because it can protect against various diseases generated by reactive oxygen species (ROS), such as neurodegenerative diseases [16], [17], cardiometabolic diseases and cancers [18], [19], [20].

Quite recently, compounds have been postulated for the inhibition of the main protease of the virus with a specific interaction in the potential treatment of infected patients [21], [22]. However, secondary metabolites are an unexplored source of chemical resources to control SARS-CoV-2 proliferation. Among natural products there are substances such as nitrogenous, terpenic and phenolic compounds, and some of them exhibit antiviral properties [23], [24]. In the polyphenol family we find flavonoids, and literature reports show the ability of some of them, such as formononetin and penduletin, to inhibit virus replication. Both are reported to be effective in vivo models against enterovirus infections that cause hand, foot and mouth disease [25]. In addition, baicalin inhibits replication of the Dengue virus in an in vitro model [26]. The increased effect of the biflavone amentoflavone compared to flavones and biflavonoid derivatives with methoxy groups [27], and the enzyme inhibiting activity of pectolinarin, herbacetin and rhoifolin [28], have been reported to inhibit SARS-CoV 3CLpro. The flavonoids hesperidin and naringin, which are abundant in citrus fruits, exhibit potential in the inhibition of the proteins responsible for the replication and propagation of SARS-CoV-2 [29]. In silico studies show that caflanone, hesperetin, myricetin and 5′-chloroquercetin can bind with high affinity to the spike protein, helicase, and protease sites on the ACE2 receptor [30], and penimethavone, a molecule mined from microorganisms, presents a high potential to modulate/inhibit the SARS-CoV-2 Mpro active site [31].

In this work, we present an in silico investigation of the antioxidant capacity and of the SARS-CoV-2 Mpro inhibitory activity of 30 polyphenolic molecules derived from Theobroma cacao: flavonoids, hydroxybenzoic acids, hydroxycinnamic acids and N-phenylpropenoyl-L-amido acids. The computational study started by evaluating some reactivity descriptors, defined in the framework of density functional theory (DFT), in order to inspect the reactive nature of these molecules. This methodology has been tested and applied in the literature by several research groups, proving to be very useful to rationalize the reactivity patterns of molecular systems.

Subsequently, the evaluation of 30 flavonoids by molecular docking showed promising results on the anti-SARS-CoV-2 activity of 5 compounds: Amentoflavone, isorhoifolin, nicotiflorin, naringin, and rutin. These compounds were subjected to more specific computational analysis to establish their potential inhibitory activity. A set of molecular dynamics simulations and free energy calculations, using MMGBSA and quantum mechanical calculations, established that all the selected compounds have binding energies very close to that of N3, a reference molecule with a high inhibitory activity against SARS-CoV-2.

We emphasize that isorhoifolin and rutin possess more negative binding energies than N3, approximately by 5 kcal mol−1. These results can be associated with the high proximity and formation of hydrogen bonds with Glu166, Cys145 and His41 residues, fundamental amino acids for the activity of Mpro. These promising results make it possible to propose these selected compounds for future experimental tests.

Computational methods

Density functional theory (DFT) calculations

Thirty polyphenolic compounds selected for their possible capability to inhibit SARS-CoV-2 protease Mpro, and N3 (reference compound) obtained from the Protein Data Bank (PDB) [32] (PDB id: 6LU7), were drawn using Discovery Studio [33] 3.1 (Accelrys, CA) and geometries were optimized using the M06–2X-D3 method [34], [35] in conjunction with the 6–31G(d,p) basis set. M06–2X-D3 is the best dispersion-corrected meta-GGA hybrid functional on the GMTKN30 database [36], and it is implemented in the Gaussian16 software program suite [37]. The water was simulated as a solvent using the SMD parametrization of the IEF-PCM.

Some DFT-based global reactivity descriptors ( Table 1), such as electronegativity (χ), global hardness (η), electrophilicity (ω), electrodonating (ω ), electroaccepting (ω ) and net electrophilicity (Δω ), were calculated to better understand the molecules reactivity.

Table 1

Equations for global reactivity indexes calculated in TAFF [38] pipeline.

	Koopmans’ theorem	Reference
Global hardness (η)	η=12ϵ_L−ϵ_H	[39], [40], [41], [42], [43], [44]
Electronegativity (χ)	χ=−12ϵ_L+ϵ_H	[40], [45], [46]
Electrophilicity (ω)	ω=μ^22η=(ϵ_L+ϵ_H)^22(ϵ_L−ϵ_H)	[47]
Electron acceptor (ω+)	ω^+=(ϵ_L+3ϵ_H)^216(ϵ_L−ϵ_H)	[47]
Electron donator (ω-)	ω^−=(3ϵ_L+ϵ_H)^216(ϵ_L−ϵ_H)	[47]
Net electrophilicity (Δω±)	∆ω^±=ω^++ω^−	[48]

Molecular docking

The 30 polyphenolic compounds and N3 were docked in the binding cavity of the main protease Mpro of SARS-CoV-2 using the AutoDock 4.0 [49] suite. The crystal structure of SARS-CoV-2 Mpro protein (2.16 Å resolution) [32] (PDB Code: 6LU7, https://www.rcsb.org/structure/6LU7) was downloaded from the PDB [50]. The SARS-CoV-2 Mpro model was modified with the Schrödinger [51] Protein Preparation Wizard [52]. Polar hydrogen atoms were added, non-polar hydrogen atoms were merged, and charges were assigned, in addition to delimiting the binding region of possible inhibitors of SARS-CoV-2 Mpro [21], [53], [54], [55], [56], [57], [58]. The binding pocket of SARS-CoV-2 Mpro was established as the center of mass between the Cys145 and His41 residues of the catalytic site. In general, the grid maps were calculated using the AutoGrid 4.0 option and the volume chosen for the grid maps was made up of 60×60×60 points, with a grid-point spacing of 0.375 Å. The author’s option was used to define the rotating bond in the ligand. In the Lamarckian genetic algorithm (LGA) dockings, an initial population of random individuals with a population size of 150, a maximum number of 2.5 × 107 energy evaluations, a maximum number of generations of 27,000, a mutation rate of 0.02 and crossover rate of 0.80 were employed. Each complex was built using the lowest docked-energy binding positions. The van der Waals interactions were computed by means of a smoothed 12–6 Lennard Jones potential, while the hydrogen bonding interactions were evaluated by a 12–10 function which incorporated a directionality term. That is, interactions which deviate from ideal hydrogen bonding geometries were progressively weighted.

The partial charges of each ligand were determined with the PM6-D3H4 semi-empirical method [59], [60] implemented in the MOPAC2016 software [61]. PM6-D3H4 introduces dispersion and hydrogen-bonded corrections to the PM6 method. 3D representations of the docking results were analyzed using the VMD molecular graphics system [62].

Ligand efficiency approach

Ligand efficiency (LE) calculations were performed using one parameter: . The parameter corresponds to the dissociation constant of a ligand/protein complex, and its value indicates the bond strength between the ligand and protein [63], [64], [65]. Low values indicate strong binding of the molecule to the protein. calculations were done using the following equations: where corresponds to binding energy (kcal mol−1) obtained from docking computations, R is the gas constant whose value is 1.987207 cal mol−1 K−1, and T is the temperature in Kelvin at standard conditions of aqueous solution at 298.15 K, neutral pH and remaining concentrations of 1 M. The ligand efficiency (LE) allows us to compare molecules according to their average binding energy [65], [66]. It is determined as the binding energy per non-hydrogen atom, as follows [63], [64], [65], [67]:where is obtained from Eq. 2 and HAC denotes the heavy atom count (i.e., number of non-hydrogen atoms) in a ligand.

Molecular dynamics simulations

MD calculations were performed for the five systems with the lowest binding energy according to the docking calculations and also for compound N3, which is our reference ligand. Each model was confined inside a periodic simulation box. The compounds were bound to SARS-CoV-2 Mpro protein [32] (PDB ID: 6LU7) in aqueous solution with an explicit solvent TIP3P water model [68] (≈20.000 water molecules). Protonation states of ionizable residues corresponding to pH 7.0 were determined by the H++ web interface that computes pK values of ionizable groups in macromolecules and adds missing hydrogen atoms according to the specified pH of the environment [69]. The N3, amentoflavone, isorhoifolin, nicotiflorin, naringin and rutin molecules were parameterized using the GAFF Force Field for organic molecules [70], [71], using the Antechamber module in AmberTools18.

The partial charges of each compound were determined by the restrained electrostatic potential (RESP) model [72] M06–2X-D3/6–31G(d,p) level. MD simulations were carried out using the modeled ff14SB [73] force field [74], [75] within the AMBER-GPU Implementations18 [76]. The simulations were carried out using a standard MD protocol: i) Minimization and structural relaxation of water molecules with 2000 steps of minimization and MD simulation with an NPT (300 K) assembly for 1000 ps using harmonic restrictions of 10 kcal mol Å−2 for protein and ligand; ii) minimization of the complete structure considering 6500 steps of conjugate gradient minimization; iii) the minimized systems were progressively heated to 300 K, with harmonic restrictions of 10 kcal mol Å−2 in the backbone protein and ligand during 0.5 ns; iv) the system was then equilibrated for 0.5 ns maintaining the restrictions and then for 5 ns without restrictions at 300 K in a canonical assembly (NVT); v) finally, the total duration of simulation was approximately 75 ns for each system. During the MD simulations, motion equations were integrated with a 2 fs time step in the NPT ensemble at a pressure of 1 atm. The SHAKE algorithm was applied to all hydrogen atoms, and the van der Waals cutoff was set to 12 Å. The temperature was maintained at 310 K, employing the Langevin thermostat method with a relaxation time of 1 ps. The Berendsen barostat was used to control the pressure at 1 atm. Long-range electrostatic forces were taken into account by means of the particle-mesh Ewald (PME) approach. Data were collected every 1 ps during the MD runs. Molecular visualization of the systems and MD trajectory analysis were carried out with the VMD software package [62].

Free energy calculation

The molecular MM-GBSA method was employed to estimate the binding free energy of the SARS-CoV-2 Mpro/ligand complexes. From calculations for a total of 75 ns of MD, the last 70 ns were extracted for analysis, and the explicit water molecules and ions were removed. The MM-GBSA analysis was performed on three subsets of each system: the protein alone, the ligand alone, and the complex (protein-ligand). For each of these subsets, the total free energy () was calculated as follows:where is the bonded (bond, angle and dihedral) and non-bonded (electrostatics and Lennard-Jones) terms; is the polar contribution of solvation energy and non-polar contribution to the solvation energy; T is the temperature; and is the conformational entropy [77]. Both and were calculated using AMBER-GPU Implementations18 [76] software with the generalized Born implicit solvent model [78], [79]. was calculated as a linear function of the solvent-accessible surface area, which was calculated with a probe radius of 1.4 Å [80]. The binding free energies of SARS-CoV-2 Mpro and ligand complexes () were calculated as differences where values are the average values over the simulation.

Non-covalent interactions

The principal clusters of the main component analyses of trajectory were analyzed with the non-covalent interaction index (NCI) [81], [82] using the NCIPLOT program to identify and map non-covalent interactions, such as hydrogen bonds, steric repulsion, and van der Waals interactions, a M06–2X-D3/6–31G(d,p). For the N3 and the antioxidants and for the protein-ligand complexes we used promolecular densities (), computed as the sum of all atomic contributions. The NCI is based on the electron density (ρ), its derivatives, and the reduced density gradient (s). The reduced density gradient is given by:

These interactions are local and manifest in real space as low-gradient isosurfaces with low densities which are interpreted and colored according to the corresponding values of sign(λ 2)ρ, where λ 2 represents the second eigenvalue of the Hessian matrix. The surfaces are colored on a blue-green-red scale according to the strength and type of interaction. Blue indicates strong attractive interactions, green indicates weak van der Waals interactions, and red indicates a strong non-bonded overlap.

Results and discussion

Global reactivity molecular descriptors

The global reactivity indices, approximated using the Koopmans’ theorem, are evaluated to predict the ligand molecules reactivity against and biological receptors. The electronegativity (χ), global hardness (η) electrophilicity (ω), electrodonating (ω ), electroaccepting (ω ) and net electrophilicity (Δω ) of antioxidants from Theobroma cacao L. and the reference molecule N3 are presented in Table 2. In general terms we can see that amentoflavone, isorhoifolin, nicotiflorin, naringin and rutin compounds have electronegativity values very close to N3′s. Now, systems with high electronegativity values tend to accept H-bonds, and the bond strength should increase with the electronegativity of the two bonded atoms, i.e., hydrogen bonding interactions could form in the active site. On the other hand, these molecules have lower hardness values than the reference molecule N3, therefore they are less averse to the arrival of electrons, and thus they are more reactive systems. These results contribute to proposals in the literature that compounds with phenol groups are more susceptible to electron donation [83]. In the case of electrophilicity values, there is an increase relative to the value of N3, indicating that these molecules are more electrophilic.

Table 2

Global reactivity descriptors for the ligands in SARS-CoV-2 Mpro complexes, calculated with the M06–2X-D3 density functional. All units measured in eV.

Compound	χ	η	ω	ω-	ω+	Δω±
N3	4.45	3.90	2.54	5.26	0.80	6.06
Amentoflavone	4.41	3.10	3.13	5.72	1.31	7.03
Naringin	4.45	3.38	2.93	5.58	1.12	6.71
Isorhoifolin	4.48	3.15	3.19	5.82	1.34	7.17
Rutin	4.26	3.13	2.89	5.42	1.15	6.57
Nicotiflorin	4.40	3.07	3.15	5.74	1.33	7.07
Kaempferol-7-O-neohesperidoside	4.20	2.98	2.96	5.44	1.23	6.68
Prunin	4.31	3.46	2.68	5.27	0.95	6.23
Kaempferol 3-O-β-D-glucoside	4.25	3.20	2.82	5.35	1.09	6.44
Isoquercetin	4.23	3.12	2.86	5.37	1.13	6.51
Luteolin 7-O-β-D-glucoside	4.34	3.18	2.96	5.53	1.19	6.73
Apigenin 7-O-glucoside	4.34	3.25	2.90	5.48	1.14	6.62
Hyperoside	4.36	3.07	3.09	5.66	1.30	6.96
Quercitrin	4.34	3.04	3.10	5.65	1.31	6.96
Isoorientin	4.36	3.12	3.04	5.61	1.25	6.87
Chlorogenic acid	4.42	3.06	3.19	5.79	1.36	7.16
Luteolin	4.28	3.18	2.88	5.41	1.13	6.55
Quercetin	4.28	2.87	3.19	5.69	1.41	7.11
Orientin	4.41	3.09	3.15	5.74	1.33	7.08
Epicatechin	3.70	3.80	1.80	4.13	0.42	4.56
(+)-Catechin	3.63	3.83	1.72	4.02	0.38	4.40
(-)-Catechin	3.68	3.75	1.81	4.12	0.43	4.56
Apigenin	4.44	3.16	3.12	5.74	1.29	7.04
Naringenin	4.24	3.54	2.54	5.10	0.86	5.97
Isorhamnetin	4.13	2.91	2.92	5.35	1.22	6.57
Isovitexin	4.38	3.19	3.00	5.60	1.21	6.81
Caffeic acid	4.35	3.07	3.08	5.64	1.28	6.93
Ferulic acid	4.33	3.06	3.05	5.60	1.27	6.88
Gallic acid	4.21	3.59	2.47	5.03	0.81	5.84
Protocatechuic acid	4.25	3.63	2.49	5.08	0.82	5.90
Coumaric acid	4.62	3.20	3.33	6.05	1.42	7.47

NCI plots at DFT geometric optimization of N3, amentoflavone, isorhoifolin, nicotiflorin, naringin and rutin at the M06–2X-D3/6–31G(d,p) level reveal different intramolecular interactions ( Fig. 1) in terms of polar groups, which can contribute to the stability in the active site of a protein through the formation of non-covalent interactions such as hydrogen bonds and ionic interactions. Note that the regions corresponding to the hydrogen bridges (blue) are tiny and surrounded by repulsive interactions (red).

Fig. 1

NCIplot of the non-covalent interaction regions with isosurface gradient (0.7 au) for A) N3, B) Amentoflavone, C) Isorhoifolin, D) Nicotiflorin, E) Naringin and F) Rutin molecules. The arrows indicate potential intramolecular interactions. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Molecular docking analysis

Mpro is a homodimeric cysteine protease and plays an important role in the SARS virus replication and transcription. When the mRNA of the virus is translated into polyproteins, Mpro is first self-cleaved to become a mature enzyme, which in turn cleaves all of the 11 remaining downstream nonstructural proteins of the polyproteins into polypeptides which are required for the replication process of the virus [84]. Consequently, Mpro is a key drug target for the inhibition of SARS-CoV-2. The binding site of Mpro is made up of subsites S1, S2, S3, S4, and S1′, which are represented based on the binding position of the substrate polyprotein [57]. The binding site is located in the groove between a His41 and Cys145 catalytic dyad. Cys145 is involved in the covalent binding of some Mpro inhibitors, specifically N3. Tung Ngo et al. reported in recent studies that Glu166 also has a prominent and important role in binding ligands to Mpro [85]. In order to investigate the possible mechanism by which selected flavonoids act, molecular docking analysis of all the flavonoids was done in the active site of Mpro. The flavonoid names and docking scores are displayed in Table 3. Five prominent ligands with the highest affinities for the active site are amentoflavone, naringin, isorhoifolin, rutin and nicotiflorin, with binding energies of −10.0, −9.0, −8.8, −8.7 and −8.5 kcal mol−1, respectively. These values are more negative than those obtained for N3 (7.5 kcal mol−1). Furthermore, it is clearly shown that glycosylated flavonoids show the highest scores, with a maximum value of 10.0 kcal mol−1, which is complemented by the study carried out by Cherrak et al. [86]. Amentoflavone showed hydrogen bond interactions with Asn142, Glu166 and Thr190. Naringin showed hydrogen bonds with His163, His164 and Phe140. Isorhoifolin showed hydrogen bonds with Ser305, Cys145 and Thr24. Rutin showed hydrogen bonds with Glu166 and Gln189. Nicotiflorin showed hydrogen bonds with Asn142, Glu166 and Cys44. The reference ligand N3 shows four hydrogen bonds with Arg188, Gln189, Cys145 and Gly143 residues. This allows us to conclude that the hydrogen bonds and the hydrophobic interactions dominate formation of these complexes.

Table 3

Molecular docking study between selected ligands and SARS-CoV-2Mpro. Intermolecular docking values, presented with their interaction energy (∆E), H-bond residues, interacting residues, and Ligand efficiency calculation for SARS-CoV-2Mpro complexes are shown.

Compound	Docking results		Ligand efficiency
	∆Ebinding (kcal mol−1)a	Residue interactionsb	Kd	LE kcal mol−1)
N3c	-7.5	ARG188, ASN142, ASP187, CYS145, GLN189, GLU166, GLY143, HIS164, HIS163, HIS41, LEU141, LEU167, MET165, MET49, PHE140, PRO168, SER144, SER305, SER46, THR25	3.18 × 10−6	0.15
Amentoflavone	-10.0	ASN142, CYS44, GLN189, GLU166, HIS41, MET165, MET49, PRO168, SER46, THR190, THR25	4.69 × 10−8	0.25
Naringin	-9.0	ASN142, CYS145, GLU166, GLY143, HIS164, HIS163, HIS41, LEU141, MET165, MET49, PHE140, SER144, SER305, SER46, THR24, THR25, THR26, THR45	2.53 × 10−7	0.21
Isorhoifolin	-8.8	ASN142, CYS145, GLN189, GLU166, GLY143, HIS163, HIS172, MET165, MET49, PHE140, SER144, SER305, SER46, THR24, THR25, THR45	3.55 × 10−7	0.21
Rutin	-8.7	CYS145, GLN189, GLU166, GLY143, GLY170, HIS163, HIS172, LEU167, MET165, MET49, PHE140, PRO168, SER144, SER305, THR25	4.20 × 10−7	0.20
Nicotiflorin	-8.5	ARG188, ASN142, CYS145, GLN189, GLU166, GLY143, HIS164, HIS163, HIS41, LEU141, LEU27, MET165, MET49, PHE140, SER144, SER305, SER46, THR25	5.89 × 10−7	0.20
Kaempferol-7-O-neohesperidoside	-8.2	ASN142, CYS145, CYS44, GLU166, GLY143, HIS164, HIS163, LEU141, MET165, PHE140, SER144, SER305, SER46, THR24, THR25, THR26, THR45	9.78 × 10−7	0.19
Prunin	-8.2	CYS145, CYS44, GLU166, HIS163, MET165, MET49, PHE140, SER144, SER46, THR25, THR26, THR45	9.78 × 10−7	0.26
Astragalin	-8.1	ARG188, ASN142, CYS44, GLN189, GLU166, HIS164, HIS163, HIS41, LEU141, MET165, MET49, PHE140, SER144, SER305, SER46, THR25, THR45	1.15 × 10−6	0.25
Isoquercetin	-8.1	ARG188, ASN142, CYS145, GLN189, GLU166, GLY143, HIS163, LEU141, MET49, PHE140, SER144, SER305, THR190, THR25	1.15 × 10−6	0.24
Luteolin 7-O-β-D-glucoside	-8.1	CYS145, GLU166, HIS163, HIS41, LEU141, MET165, MET49, PHE140, SER144, SER305, SER46, THR24, THR25, THR45	1.15 × 10−6	0.25
Apigenin 7-O-glucoside	-7.9	CYS145, GLU166, HIS164, HIS41, LEU141, MET165, MET49, PHE140, SER144, SER305, SER46, THR24, THR25, THR26, THR45	1.62 × 10−6	0.25
Hyperoside	-7.7	ASN142, CYS145, GLN189, GLU166, GLY143, HIS163, HIS172, LEU141, MET49, PHE140, SER144, SER46, THR25, THR26	2.27 × 10−6	0.23
Quercitrin	-7.7	ASN142, CYS145, GLN189, GLU166, GLY143, HIS164, HIS163, HIS172, LEU141, MET165, MET49, PHE140, SER144, SER305, THR25	2.27 × 10−6	0.24
Isoorientin	-7.6	ARG188, ASN142, CYS145, GLN189, GLU166, GLY143, HIS41, MET165, MET49, SER46, THR24, THR25	2.69 × 10−6	0.23
Chlorogenic acid	-7.4	ASN142, CYS145, CYS44, GLU166, HIS163, LEU141, MET165, MET49, PHE140, SER144, THR25, THR45	3.77 × 10−6	0.29
Luteolin	-7.4	ARG188, GLN189, GLU166, HIS164, HIS163, LEU141, MET165, MET49, PHE140, SER144	3.77 × 10−6	0.35
Quercetin	-7.3	CYS145, GLU166, HIS164, HIS163, HIS41, LEU141, MET165, MET49, PHE140, SER144, SER305, THR25	4.46 × 10−6	0.33
Orientin	-7.3	ALA191, ASN142, GLN189, GLU166, HIS164, HIS41, MET165, MET49, PRO168, THR190	4.46 × 10−6	0.22
Epicatechin	-7.2	ARG188, GLN189, GLU166, HIS163, LEU141, MET165, MET49, PHE140, SER144	5.28 × 10−6	0.34
(+)-Catechin	-7.1	ARG188, CYS145, GLN189, GLU166, HIS164, HIS41, LEU141, MET165, MET49, PHE140, SER305	6.26 × 10−6	0.33
(-)-Catechin	-7.1	ASN142, GLN189, GLU166, HIS164, HIS163, HIS172, HIS41, MET165, MET49, PHE140, SER305	6.26 × 10−6	0.33
Apigenin	-7.1	ARG188, GLN189, GLU166, HIS163, LEU141, MET165, MET49, PHE140, SER144	6.26 × 10−6	0.35
Naringenin	-6.9	ARG188, GLN189, GLU166, HIS163, LEU141, MET165, MET49, PHE140, SER144	8.77 × 10−6	0.34
Isorhamnetin	-6.9	CYS145, GLU166, HIS163, HIS41, LEU141, MET165, MET49, PHE140, SER144, SER305, THR25	8.77 × 10−6	0.30
Isovitexin	-6.9	ARG188, GLN189, GLU166, HIS164, HIS41, MET165, MET49, SER46, THR24, THR25, THR45	8.77 × 10−6	0.22
Caffeic acid	-5.4	CYS145, GLU166, HIS164, HIS163, LEU141, MET165, MET49, PHE140, SER144, SER305	1.10 × 10−4	0.41
Ferulic acid	-5.3	CYS145, GLU166, HIS164, HIS163, HIS41, LEU141, MET165, PHE140, SER305	1.30 × 10−4	0.37
Gallic acid	-5.2	ASN518, GLY170, GLY306, LEU586, LYS137, PHE307, SER305, THR169, VAL171	1.54 × 10−4	0.43
Protocatechuic acid	-5.2	ARG308, GLY138, GLY170, GLY306, LEU586, LYS137, PHE307, THR169, VAL171	1.54 × 10−4	0.47
Coumaric acid	-5.2	ASN142, GLU166, HIS163, HIS41, LEU141, MET165, MET49, PHE140, SER144	1.54 × 10−4	0.43

*Residues His41 and Cys145 of the catalytic site are highlighted with bold font.

In each site, the energy was calculated to see which bound more tightly to the ligand.

A cutoff distance of 3 Å was chosen to consider that a hydrogen bond was established.

N3 is our reference ligand and the structure of its complex with Mpro was obtained from PDB (id: 6LU7).

Ligand efficiency analysis

The parameters K and LE were used to compare the affinity of the molecules bound to SARS-CoV-2 Mpro in this study. K corresponds to the dissociation constant of a ligand-protease complex. Therefore, very low values indicate that the compound binds tightly to the protein. LE represents the average binding energy per non-hydrogen atom, where tolerable values of LE for inhibitor candidates are LE > 0.3 kcal mol−1 (see Table 2). The five best ligands, amentoflavone, isorhoifolin, nicotiflorin, naringin, rutin and the reference molecule N3, as obtained from the docking, exhibit low K values, which means that these complexes are the most stable in the series. Although the LE values are close to 0.3 kcal mol−1, this does not affect binding efficiency to the active site of the protease. The results are consistent with those obtained in the molecular docking in which these complexes were the most stable according to their ∆E values. Therefore, these molecules are excellent prospects for testing as SARS-CoV-2 Mpro inhibitors.

Molecular dynamics simulation and MM/GBSA analysis

Molecular dynamics simulations were performed for 75 ns to analyze the steady nature and conformational stability of SARS-CoV-2 Mpro bound to N3, amentoflavone, isorhoifolin, nicotiflorin, naringin and rutin.

The RMSD results along the MD trajectories show high stability for the SARS-CoV-2 Mpro dimer, as illustrated in Fig. 2. After a molecular dynamics simulation of 75 ns the N3, amentoflavone, isorhoifolin, nicotiflorin and rutin molecules remain within parameters that agree wity the system being in equilibrium. Therefore, no complex suffered structural destabilization during the simulation, except for the naringin complex with the SARS-CoV-2 Mpro. Instability is observed from 42ns where the RMSD values oscillate by more than 3.0 Å (see Fig. 2-E). Deviations with a maximum difference of 3.0 Å of RMSD [87] indicate that the system is in equilibrium. Also, the RMSD curves for N3 and amentoflavone are remarkably more stable than for isorhoifolin, nicotiflorin and rutin, with RMSD values differing by about 1.0 Å between the above-mentioned complexes.

Fig. 2

Root Mean Square Deviation (RMSD) as a function of simulation times for the complexes formed between SARS-CoV-2 Mpro and A) N3, B) Amentoflavone, C) Isorhoifolin, D) Nicotiflorin, E) Naringin and F) Rutin.

The Radius of Gyration (RGyr) study complements the RMSD analysis. RGyr was computed for the same runs. The results shown in the boxplot (see Fig. 3) illustrate the distribution of the data within RGyr as a function of time. The trajectories show that the RGyr of N3 and amentoflavone have compactly distributed interquartile range values, showing uniform fluctuations of RGyr as a function of the molecular trajectory. Their mean values are 4.97 and 4.74 Å. In the case of isorhoifolin, rutin, nicotiflorin and naringin, the Interquartile Range values are widely distributed, meaning that their RGyr values fluctuate more depending on the molecular trajectory. Their mean values are 5.31, 4.51, 4.25 and 4.84 Å. Although the fluctuation of the interquartile range of RGyr for these cases is higher, these never exceed the N3 and amentoflavone values byas much as 0.6 Å, meaning that the ligands remain quite stable during the 75 ns simulation for RGyr, indicating good thermodynamic stability of their Mpro binding, with only minor conformational changes.

Fig. 3

Box-plot displaying radii of gyration for the SARS-CoV-2 Mpro in complex with N3, amentoflavone, isorhoifolin, nicotiflorin, naringin and rutin.

After verifying that SARS-CoV-2 Mpro maintains its folding when bound to the ligands, it is desirable to perform a quantitative analysis of the structural fluctuation at the residue level using the Cα atom of each amino acid. To compare in detail the differences between the complexes a suitable measure is Normal Mode Analysis (NMA) and Root Mean Square Fluctuation (RMSF) values, which indicate the fluctuation of the structure with respect to the average conformation in the 75 ns simulation. In this way, the conformational stability of the structure is quantified. When overlaying the RMSF profiles of SARS-CoV-2 Mpro with different ligands (top Fig. 4), the Mpro complex present with naringin showed the greatest deviation of its residues. In contrast, the other complexes show less fluctuation, which could mean greater stability of the backbone in the Mpro complexes with amentoflavone, isorhoifolin, nicotiflorin and rutin. The most notable differences are found in the Tyr54, Glu55, Ala193, Ala194 and Gly195 residues, which show higher RMSF values. In the active site of the protease, the fluctuation values of the main residues [21] (His41, His163, His164, Phe140 and Cys145) were similar among the complexes with the five selected molecules. In the case of NMA profiles (bottom Fig. 4), they showed similar trends for Mpro bound to isorhoifolin, nicotiflorin and rutin, and different trends with amentoflavone and naringin. The amentoflavone-Mpro and naringin-Mpro complexes showed a comparatively larger fluctuation of residues near Asn51 and the terminal residues of Mpro. This larger displacement in their normal modes did not affect the flexibility of the protein throughout the simulations. These results indicate that isorhoifolin and rutin are very likely to have the same capability to inhibit the SARS-CoV-2 Mpro as N3, but reversibly.

Fig. 4

NMA and RMSF of the α-carbons. A principal component analysis was carried out using the 75 ns trajectories, and the main normal mode of movement was obtained. The displacement was plotted for each residue of SARS-CoV-2 Mpro in complex with N3, amentoflavone, isorhoifolin, nicotiflorin, naringin and rutin. Pocket site residues are distinguished in gray boxes. His41 and Cys145 residues of the catalytic site are highlighted with bold font.

The analyses of intermolecular interactions performed during the molecular simulation show that N3, amentoflavone, isorhoifolin, nicotiflorin, naringin and rutin preserve hydrogen bonds with active site residues of the SARS-CoV-2 Mpro. However, the number of hydrogen bonds formed was different for each ligand ( Fig. 5, Fig. 6). N3 formed 3 hydrogen bonds between the residues Thr190, Gln192 and Glu166, highlighting the participation of the residues Gln189, Met49, His41 and Cys145. Amentoflavone formed 4 hydrogen bonds with Cys44, Glu166, His41 and Thr190, highlighting the participation of Asn142, Ser46 and Met49. In the case of isorhoifolin, 3 hydrogen bonds between es Asp187, Glu166 and Leu141 were determined, highlighting the participation of Gln189 and Phe140. Nicotiflorin formed 4 hydrogen bonds with Asp187, Glu166, Asn142 and Gln189, highlighting the participation of Arg188, His164, Ser46, Cys145 and His41. Naringin formed 2 hydrogen bonds with Glu166 and Asn142, highlighting the participation of ues Glu47, Gln192 and Gln189. Three hydrogen bonds are formed between rutin and Glu166, Thr25 and Asp187, highlighting the participation of Cys145, Gln189, His163 and His41. Finally, these residues (see Fig. 5, Fig. 6) are consistent with previous theoretical-experimental studies carried out by Dai et al. [55], in which they detail the interactions that some of the synthesized compounds have with the active site of the Mpro.

Fig. 5

Frequency of the appearance of residues at a distance of 3.0 Å or less from ligands A) N3, B) Amentoflavone, C) Isorhoifolin, D) Nicotiflorin, E) Naringin and F) Rutin calculated using MD procedures.

Fig. 6

Fraction of intermolecular hydrogen bonds for SARS-CoV-2 Mpro interacting with A) N3, B) Amentoflavone, C) Isorhoifolin, D) Nicotiflorin, E) Naringin and F) Rutin. The graph bar shows the most common hydrogen bonds formed between the residues in the binding pocket and the inhibitors. Values obtained from CPPTRAJ script in AMBER.

The binding mode of the molecules in the active site of the SARS-CoV-2 Mpro at equilibrium is displayed in Fig. 7. The NCIplot analysis labeled all the hydrogen-bonding and van der Waals interactions in total agreement with the classical molecular dynamics simulations, providing a qualitative confirmation of these interactions, using a topological and visual analysis of a scalar field related to the electron density and reduced density gradient (Fig. 7-III). These results suggest that the improved character of the inhibitors is due to direct interactions with the protease.

Fig. 7

Schematic representations of principal component analysis of the respective production runs for A) Amentoflavone, B) Isorhoifolin, C) Nicotiflorin, D) N3, E) Naringin and F) Rutin bound to SARS-CoV-2 Mpro. (I) Representative amino acid residues surrounding ligands in the binding pocket of Mpro. (II) Two-dimensional interaction map of ligands and Mpro. (III) NCIPLOT isosurface gradient (0.5 au) of ligands in the structure of Mpro. The arrows indicate potential interactions between amino acid residues and ligands.

The binding free energy (MM-GBSA) was calculated after MD simulation. The last 70 ns were collected for all complexes, and the results are shown in Table 4. Isorhoifolin and rutin showed better binding free energies with values of −43.91 and −44.91 kcal·mol−1, rutin being the molecule with the lowest binding free energy, while amentoflavone, nicotiflorin and naringin showed less negative binding free energies with values of −33.77, −34.43 and −19.10 kcal mol−1, respectively. Furthermore, N3 showed a binding free energy value of −38.87 kcal mol−1 with Mpro of SARS-CoV-2, this being the third best value compared to isorhoifolin and rutin, with energy differences of 5.04 and 6.04 kcal mol−1. Although 3N is a synthetic molecule with Mpro inhibitory activity for SARS-CoV-2, the binding free energy results showed two possible candidates, isorhoifolin and rutin, that are also likely to exhibit Mpro inhibitory activity. Experiments performed by Huynh and col. and Rahman and col. mention that rutin can be used as a potential Mpro inhibitor of COVID-19 [88], [89], based on its well-documented significant antibacterial and antiviral properties [28], [90], [91].

Table 4

Predicted binding free energies (kcal mol−1) and individual energy terms calculated from molecular dynamics simulation following the MM-GBSA protocol for Mpro complexes.

Compunds	Calculated free energy of decomposition (kcal mol−1)
	ΔGbinding	ΔEvdW	ΔEelect	ΔGgas	ΔGsolv
N3	-38.87 ± 0.04	-54.57 ± 0.05	-17.74 ± 0.06	-72.32 ± 0.08	33.44 ± 0.05
Amentoflavone	-33.77 ± 0.09	-41.23 ± 0.08	-25.78 ± 0.11	-67.01 ± 0.15	33.24 ± 0.08
Isorhoifolin	-43.91 ± 0.09	-45.54 ± 0.08	-62.18 ± 0.16	-107.72 ± 0.17	63.81 ± 0.13
Nicotiflorin	-34.43 ± 0.10	-43.99 ± 0.07	-41.23 ± 0.22	-85.23 ± 0.22	50.80 ± 0.14
Naringin	-19.10 ± 0.12	-30.97 ± 0.13	-14.64 ± 0.19	-45.62 ± 0.25	26.51 ± 0.15
Rutin	-44.91 ± 0.08	-55.18 ± 0.06	-42.78 ± 0.19	-97.97 ± 0.18	53.06 ± 0.11

N3, rutin and isorhoifolin as antiviral molecules

Consistent with the results obtained in this work, rutin has been reported to show antimicrobial activity [92] and, through in silico studies, possible inhibitory activity of several proteins that are essential for SARS-CoV-2 to complete its viral cycle, [86], [88], [89], [93], [94]. However, its antiviral spectrom is broader and it is being tested experimentally as an antiviral agent against retroviruses, orthomyxoviruses, herpes viruses, hepatitis B and C viruses, and the H1N1 influenza virus. In addition to showing antiviral activity, some studies demonstrate its activity associated with anti-inflammatory processes, which is very important in patients with COVID-19 [95]. On the other hand, no antiviral activity of isorhoifolin has been shown, but anti-acaricidal activity has been reported, with an LC50 = 0.65 mg ml−1 [96], and potential inhibition of β-amyloid protein aggregation. These compounds interact with the Glu166 residue, reaching an occupancy greater than 0.6. This interaction may be relevant for its antiviral activity since Glu166 has been shown to stabilize the Mpro dimer and its mutation drastically reduces its activity [85]. Furthermore, these ligands could form hydrogen bonds with the catalytic residues Cys145 and His41, an interaction that could explain the very favorable free energy of binding of these compounds. In this sense, rutin and isorhoifolin form interactions similar to those described for the reversible stage of N3 binding [97] ( Fig. 8). As we have reported in this research, amentoflavone, isorhoifolin, nicotiflorin, naringin and rutin, all present in Criollo, Trinitario and Forastero cocoa beans [9] ( Table 5), can interact with various residues in the active site of Mpro, including catalytic residues, showing excellent potential to inhibit the replication of SARS-CoV-2. Therefore, this research is a starting point for in vitro studies and the development of pharmacological therapies against this virus, using bioproducts developed from cocoa rich in bioflavonoids such as rutin and isorhoifolin.

Fig. 8

Representative amino acid residues surrounding ligands in the binding pocket of Mpro for A) Isorhoifolin, B) Rutin and C) N3.

Table 5

Concentrations of compounds that showed the highest inhibitory activity of Mpro protease in varieties of Theobroma cacao.

No.	Compound	Class	Polyphenols (mg/kg)
			Criollo	Trinitario	Forastero
1	Amentoflavone	Flavones	n.d.	0.24	5.8
2	Naringin	Flavanones	n.d.	n.d.	n.d.
3	Isorhoifolin	Flavones	37.5	95	44.7
4	Rutin	Flavonols	50	20	10
5	Nicotiflorin	Flavonols	12.1	7.4	8.5
6	Kaempferol-7-O-neohesperidoside	Flavonols	9.9	2.85	3.9
7	Prunin	Flavanones	33.8	54.65	31.6
8	Astragalin	Flavonols	0.5	0.5	0.4
9	Isoquercetin	Flavonols	n.d.	n.d.	n.d.
10	Luteolin 7-O-β-D-glucoside	Flavones	1.2	0.55	1.0

Conclusions

Theobroma cacao contains polyphenolic compounds, secondary metabolites that are widely available in natural medicines. The identification of these polyphenolic compounds is a remarkable step towards the discovery of drugs targeting SARS-Cov-2 Mpro, and may help in the development of new drugs. Results suggest that the flavonoids isorhoifolin and rutin should display strong inhibitory activities against SARS-CoV-2 Mpro. Isorhoifolin and rutin seem to bind more strongly than N3 co-crystallized inhibitor. The estimated free energies of binding of these two flavonoids are −43.91 and −44.91 kcal mol−1, respectively. The binding of the flavonoids at the protease active site is structure dependent, proffering the spontaneous and energetically favored production of the protein ligand complex. The interactions of flavonoids with active site residues mainly occur through the hydroxyl groups in their structure. These results represent an important step in the exploration of natural products from Theobroma cacao for structure-based design of anti-SARS-CoV-2 drugs. Furthermore, they encourage further in vitro research and also boost the traditional preventive use of isorhoifolin and rutin. We anticipate that the insights obtained in the present study may prove valuable for researching and developing novel anti-SARS-CoV-2 coronavirus therapeutic agents in the future.

CRediT authorship contribution statement

Osvaldo Yañez, Olimpo García-Beltrán, Carlos Areche, William Tiznado and Angelica Sandoval-Aldana contributed to the conception and design of the study; Osvaldo Yañez, Alejandro Vasquez-Espinal, Edison Osorio, Manuel Isaías Osorio and Jessica Bravo preformed the theoretical calculations; Osvaldo Yañez, José M. Pérez-Donoso, Fernando González-Nilo, Manuel Isaías Osorio and Angelica Sandoval-Aldana organized the database; Osvaldo Yañez, William Tiznado, Carlos Areche and Olimpo García-Beltrán wrote the first draft of the manuscript; Osvaldo Yañez, Olimpo García-Beltrán and Maria João Matos wrote sections of the manuscript. All authors contributed to manuscript revision, read and approved the submitted version.

Declaration of Competing Interest

No conflict of interest.