DFT, molecular docking and molecular dynamics simulation studies on some newly introduced natural products for their potential use against SARS-CoV-2.

Throughout the history, natural products always give new paths to develop new drugs. As with many other diseases, natural compounds can be helpful in the treatment of COVID-19. SARS-CoV-2 main protease enzyme has an important role in viral replication and transcription. Therefore, inhibiting this enzyme may be helpful in the treatment of COVID-19. In this study, it is aimed to investigate eight natural compounds which have recently entered the literature, computationally for their potential use against SARS-CoV-2. For this purpose, first, density functional theory (DFT) calculations were performed on the investigated compounds, and energy minimizations, geometry optimizations, vibrational analyses, molecular electrostatic potential map calculations were carried out. After DFT calculations, geometry optimized structures were subjected to molecular docking calculations with the use of SARS-CoV-2 main protease (pdb id: 5r80) and top-scoring ligand-receptor complexes were obtained. In the next part of the study, molecular dynamics (MD) simulations were performed on the top-scoring ligand-receptor complexes to investigate the stability of the ligand-receptor complexes and the interactions between ligands and receptor in more detail. Additionally, in this part of the study, binding free energies are calculated with the use of molecular mechanics with Poisson-Boltzmann surface area (MM-PBSA) method. Results showed that, all ligand-receptor complexes remain stable during the MD simulations and most of the investigated compounds but especially two of them showed considerably high binding affinity to SARS-CoV-2 main protease. Finally, in the study, ADME (adsorption, desorption, metabolism, excretion) predictions and drug-likeness analyses were performed on the investigated compounds.

Introduction

Since its first appearance in December 2019, SARS-CoV-2 continues to spread rapidly all over the world. Although the number of deaths due to COVID-19 is dramatically increasing, an effective drug treatment other than vaccines has not been developed yet. The methods currently studied for the treatment of COVID-19 can be classified in three groups: (1) Vaccine development studies, (2) New drug development studies, (3) Drug repurposing studies. In new drug development and drug repurposing studies SARS-CoV-2 main protease is a commonly selected target because of its important role in viral replication and transcription. Thus, inhibiting this enzyme may be a potential treatment of COVID-19. The extant literature contains several studies on SARS-CoV-2 main protease and its potential inhibitors [1], [2], [3], [4], [5], [6], [7], [8]. Throughout the history, natural compounds have been used in the treatment of various diseases, and today they continue to be important compounds in the development of many drugs.

In this study, eight natural compounds which have recently been added to the literature by Arriffin et al. [9], Zhang et al. [10], Liu et al. [11] and Yao et al. [12] have been investigated computationally for their potential to inhibit the SARS-CoV-2 main protease. In the study, DFT calculations, molecular docking calculations, molecular dynamics simulations, binding free energy calculations and finally drug-likeness and ADME analyses were performed on the investigated compounds. Result showed that investigated compounds, especially compounds 1 and 7, showed considerably high binding affinity to SARS-CoV-2 main protease and may be promising structures for the treatment of COVID-19. Chemical structures of the investigated compounds are given in Fig. 1 [9], [10], [11], [12].

Fig. 1

Chemical structures of the investigated compounds.

Materials and methods

Density functional theory (DFT) calculations

DFT calculations were performed with the use of Becke three-parameter hybrid functional combined with Lee-Yang-Parr correlation functional (B3LYP) and 6-311+G(2d,p) basis set. Calculations were performed with the use of IEFPCM solvation model and water was selected as solvent. In this part of the study, geometry optimizations, frequency analyses, molecular electrostatic potential (MEP) map calculations were carried out. Prior to geometry optimizations, a conformational search was performed for each structure. In this part, a frequency analysis was also carried out for each structure to confirm that all geometry optimized structures correspond to a real minimum. All calculations in this part were carried out with the use of Gaussian 09 Rev.D01 [13], GaussView 5 [14] and VeraChem Vconf [15] program packages, and Discovery Studio Visualizer [16] was used in the representation of the results.

Molecular docking and molecular dynamics (MD) simulation studies

Geometry optimized structures of the investigated compounds which were obtained from DFT calculations were used for molecular docking calculations. The structure of SARS-CoV-2 main protease was downloaded from RCSB PDB Data Bank (pdb id: 5r80) [17]. The docking studies were performed with the use of AutoDock Tools [18] and AutoDock Vina [19] program packages, and Discovery Studio Visualizer [16] was used in the representation of the docking results. Prior to molecular docking calculations, water molecules and the ligand in the enzyme structure were removed, hydrogens and Gasteiger charges were added. Calculations were performed with the use of Lamarckian genetic algorithm in a 20 × 20 × 20 Å3 grid box and grid spacing was taken as 1 Å. Nine docking poses were obtained for each molecular docking calculation. After molecular docking calculations, top-scoring ligand-receptor complexes were subjected to 30 ns all-atom MD simulations to investigate the ligand-receptor interactions in more detail and to determine the binding free energies accurately. In MD simulations GROMACS 5.1.4 [20] program package was used. Topology of the enzyme was prepared with the use of AMBER force field [21] and TIP3P water model. Ligand topologies were obtained from Acpype Server [22]. After neutralization of the system, energy minimization was performed for each complex by employing steepest descent minimization algorithm. After performing 200 ps NVT and NPT ensemble equilibrations, MD simulations were performed in a dodecahedron simulation box for 30 ns at 1 bar and 300 K reference pressure and temperature. After performing MD simulations, binding free energies were obtained from MM-PBSA calculations with the use of g_mmpbsa tool [23,24]. MM-PBSA calculations were performed for the last 20 ns of the MD simulations to determine the binding energies more accurately.

Drug-likeness analyses and ADME studies

Drug-likeness analyses and ADME studies are important to analyze the pharmacokinetics of a given molecule which could be used as a drug. Drug-likeness analyses and ADME predictions were performed via SwissADME web server [25]. In this part of the study, physicochemical properties, lipophilicity, water solubility and various pharmacokinetics parameters were predicted. Additionally, drug-likeness analyses were performed with the use of Lipinski [26], Ghose [27], Veber [28], Egan [29] and Muegge [30] filters, and medicinal chemistry analyses were also carried out.

Results and discussion

Optimized geometries of the investigated compounds obtained from DFT calculations are given in Fig. 2 . A frequency analysis was performed for each optimized structure to confirm that the optimized geometry corresponds to a global minimum. Results showed that all the geometry optimized structures correspond to a global minimum. Calculated geometric parameters (Tables S1–S8) and vibrational spectra (Figs. S1–S8) of the investigated compounds are given in Supplementary File.

Fig. 2

Optimized geometric structures of the compounds 1-8.

Molecular electrostatic potential (MEP) maps are very useful tools to obtain information about the electron rich and electron deficient parts of a given molecule. MEP maps of the investigated compounds obtained from DFT calculations with the use of 6-311+G(2d,p) basis set and isovalue of 0.0004 are given in Fig. 3 . It was observed that in most cases, the negative charges are located on the oxygen atoms, and almost in all cases the most positive regions of the molecules are the regions where the hydrogens bonded to oxygen atoms are located. These negative and positive centers have taken part in forming non-bonded interactions (especially hydrogen bonds) in the ligand-receptor complexes throughout molecular docking and molecular dynamics simulations.

Fig. 3

MEP maps of the compounds 1-8.

Molecular docking and molecular dynamics (MD) simulations

In molecular docking calculations nine top-scoring docked poses have been obtained for each compound. Highest docking scores for each compound are given in Table 1 . Some known drugs currently under investigation for their potential use for the treatment of COVID-19 were also subjected to molecular docking to make a comparison. 3D structures of remdesivir, favipiravir and hydroxychloroquine were obtained from PubChem Database [31] and the structure of lopinavir was obtained from Zinc Database [32]. Results showed that investigated compounds have docking scores comparable to the reference drugs.

Table 1

AutoDock Vina binding scores of the investigated compounds.

Compound	AutoDock Vina Score(kcal/mol)
1	-8.6
2	-7.5
3	-7.5
4	-7.7
5	-6.2
6	-6.6
7	-7.5
8	-7.6
Remdesivir	-7.9
Lopinavir	-8.2
Favipiravir	-5.1
Hydroxychloroquine	-6.3

After molecular docking calculations, each top-scoring ligand-receptor complexes were subjected to 30 ns MD simulations to investigate the stability of ligand-receptor complexes and to investigate the ligand-receptor interactions in more detail. Additionally, in this part, calculations with the use of MM-PBSA method were carried out to determine the binding free energies accurately. Binding positions of the investigated natural compounds and ligand-receptor interactions at the end of each 30 ns MD simulation are given in Figs. 4 , and 2D representations of these interactions are given in Fig. 5 .

Fig. 4

Binding positions and ligand-receptor interactions at the end of 30 ns MD simulations.

Fig. 5

2D representations of ligand-receptor interactions at the end of 30 ns MD simulations.

Results showed that all ligands bind to the active site of the enzyme and remain there throughout the MD simulations. It was observed that in all complexes, HIS41 and CYS145 amino acids of the SARS-CoV-2 main protease take part in the ligand-receptor interactions and in most complexes, ligands interacted with ASN142 and MET165 amino acids of the enzyme. Additionally, in four of the complexes, MET49, in three of the complexes, PRO168, GLU166, LEU167 and GLN189, and in two of the complexes, THR26, LEU27, CYS44, GLY143, ASP187 and ARG188 amino acids participated in the interactions between ligands and receptor. On the other hand, in only one complex, LEU14, HIS163, SER46, THR45, VAL186 and HIS164 amino acids took part in the interactions. Results showed that, mainly hydrogen bonds and some other type of interactions (i.e., pi-sigma, pi-pi T-shaped, alkyl, pi-alkyl, pi-pi stacked, pi-sulfur, and amide-pi stacked) take part in the stabilization of the complexes.

Molecular dynamics simulations were also performed for reference drugs, hydroxychloroquine, lopinavir, and remdesivir and favipiravir. Interactions between reference drugs and SARS-CoV-2 main protease at the end of 30 ns MD simulations are given in Fig. 6 . MD simulations showed that favipiravir could not hold its position in the binding pocket of the enzyme and detached from the enzyme.

Fig. 6

Interactions between reference drugs and SARS-CoV-2 main protease at the end of 30 ns MD simulation.

The stability of the complex during MD simulation was evaluated with the use of root mean square deviation (RMSD) and radius of gyration. RMSDs of main protease compared to its initial position were given in Fig. 7 . Average RMSDs and standard deviations were found to be 0.17 ± 0.02, 0.22 ± 0.05, 0.17 ± 0.03, 0.20 ± 0.03, 0.24 ± 0.04, 0.24 ± 0.05, 0.16 ± 0.02 and 0.19 ± 0.03 nm for the complexes of compounds 1-8, respectively. Results showed that in all ligand-receptor complexes (especially in the complexes of compounds 1, 3, 4, 7 and 8), SARS-CoV-2 main protease remained stable throughout the MD simulation. Additionally, RMSD of ligand compared to the position of the enzyme was monitored for each complex to investigate how well the binding pose was preserved during the MD simulation (Fig. 8 ). Results showed that in the complexes of compounds 1, 3 and 7, the position of the ligands preserved during the entire simulation. In the complexes of compounds 5 and 8, ligands quickly reached its equilibrium position and remained in this position during almost the entire simulation. In the complexes of compounds 2 and 4, no considerable change in the position of the ligand was observed after the first 15 ns.

Fig. 7

RMSD of SARS-CoV-2 main protease compared to its initial position.

Fig. 8

RMSD of ligand compared to SARS-CoV-2 main protease.

Radius of gyration (RG) is another parameter which is used to investigate the stability of the protein. RG of the enzyme was monitored during the MD simulation (Fig. 9 ). Average RGs and standard deviations of the enzyme in the complexes of compounds 1-8 were found to be 2.21 ± 0.01, 2.23 ± 0.02, 2.21 ± 0.01, 2.22 ± 0.01, 2.24 ± 0.02, 2.23 ± 0.01, 2.21 ± 0.02 and 2.23 ± 0.01 nm, respectively. Results showed that there is not a considerable change in the RG of enzyme, and it remained stable during the entire simulation.

Fig. 9

Changes in radius of gyration of SARS-CoV-2 main protease during the MD simulations.

RMSDs of reference drugs compared to the position of SARS-CoV-2 main protease were monitored throughout the MD simulation (Fig. 10 ). Results showed that hydroxychloroquine, lopinavir and remdesivir remained in the binding pocket throughout the entire simulations but favipiravir could not hold its position in the binding pocket. It was observed that favipiravir was detached from the enzyme structure at 10th ns and reattached at 13th ns, but it completely moved away from the enzyme at 25th ns.

Fig. 10

RMSDs of reference drugs compared to the position of SARS-CoV-2 main protease.

In order to make a comparison, RMSDs and RGs of SARS-CoV-2 main protease were also monitored for the complexes of reference drugs except favipiravir (Fig. 11 ). Average RMSDs and standard deviations for hydroxychloroquine, lopinavir and remdesivir were found to be 0.18 ± 0.02, 0.20 ± 0.03 and 0.22 ± 0.07, respectively. Additionally, average RGs and standard deviations for hydroxychloroquine, lopinavir and remdesivir were found to be 2.22 ± 0.01, 2.22 ± 0.02 and 2.23 ± 0.02, respectively. Since favipiravir could not form a stable complex with SARS-CoV-2 main protease, no further investigations were performed on the SARS-CoV-2 main protease – favipiravir complex.

Fig. 11

RMSDs and RGs of enzyme obtained for the complexes of reference drugs.

After performing MD simulations, binding free energies were obtained from molecular mechanics with Poisson-Boltzmann surface area (MM-PBSA) calculations with the use of g_mmpbsa tool [23,24]. MM-PBSA calculations were performed for the last 20 ns of the MD simulations to increase the accuracy. Average binding free energies are given in Fig. 12 . Results showed that lopinavir has the highest binding affinity to SARS-CoV-2 main protease but most of the compounds, have higher binding affinity than remdesivir, and 7, 1 and 3 have higher binding affinity than hydroxychloroquine.

Fig. 12

Average binding free energies obtained from MM-PBSA calculations (Hyd: Hydroxychloroquine, Lop: Lopinavir, Rem: Remdesivir).

Drug-likeness analyses and ADME predictions for investigated natural compounds were performed via SwissADME web server [25]. Results are given in Tables 2 and 3 . Drug-likeness analyses and ADME predictions were also performed for reference drugs, hydroxychloroquine, lopinavir, remdesivir and favipiravir to make a comparison. Results showed that none of the investigated natural compounds have more than ten hydrogen bond acceptors, in fact the highest number of hydrogen bond acceptors was obtained for compound 7 with the value of seven. On the other hand, it was observed that only compound 7 has more than five hydrogen bond donors but this is the only violation of Lipinski rule for compound 7 and does not prevent it from obeying the Lipinski rule. Molecular weights of the investigated natural compounds are in the range of 266.3-372.4 g/mol. According to Veber rule and Muegge rule, rotatable bonds in the molecule should not be more than 10 and 15, respectively. All the molecules investigated in the study have a maximum of seven rotatable bonds and obey the Veber and Muegge rules. Fraction of sp3 carbon atoms is used to characterize the aliphatic degree of a molecule and to predict solubility. It is proposed that increasing saturation increases the clinical success rate of a given molecule [33,34]. In the study, for most of the investigated compounds, fraction C sp3 was found to be higher than 0.25 and three of the molecules (compounds 5, 6 and 7) have much more higher degrees of saturation. In the study, molar refractivity of the investigated compounds was predicted to be in the range of 72.51-100.72, and it was observed that all the molecules obey Ghose rule. Topological polar surface area (TPSA) is another parameter to predict the drug-likeness of a given molecule [35]. According to Veber rule, polar surface area should not be greater than 140 Å2. In the study, it was observed that the polar surface areas of the investigated compounds are not greater than 97.99 Å2. According to Ghose rule, partition coefficient for a given molecule should be in the range of -0.4 to 5.6. In the study, the average partition coefficient of the molecules which is the average of iLOGP [36], XLOGP3 [37], WLOGP [38], MLOGP [26,39,40] and SILICOS-IT [41] were predicted to be in the range of 0.32-3.22. Water solubilities of the investigated compounds (ESOL, ALI and SILICOS-IT [41], [42], [43]) varied from moderately soluble to very soluble.

Table 2

Results of ADME analyses of the investigated compounds.

Compound	1	2	3	4	5	6	7	8	Fav	Lop	Rem	Hyd
Physicochemical Properties
Number of H-Bond Acceptors	5	5	5	6	4	5	7	4	4	5	12	3
Number of H-Bond Donors	2	2	2	1	3	4	6	3	2	4	4	2
Molecular Weight (g/mol)	338.4	354.4	300.3	314.3	266.3	284.4	372.4	284.3	157.1	628.8	602.5	335.8
Number of heavy atoms	25	26	22	23	19	20	27	21	11	46	42	23
Number of aromatic heavy atoms	12	12	12	12	0	0	12	14	6	18	15	10
Fraction C sp3	0.25	0.29	0.29	0.24	0.67	0.80	0.38	0.18	0.00	0.43	0.48	0.50
Number of rotatable bonds	1	4	1	1	3	3	7	1	1	17	14	9
Molar refractivity	92.81	98.90	80.16	80.58	72.51	74.14	100.72	83.95	32.91	187.92	150.43	98.57
TPSA (Å2)	75.99	75.99	68.15	85.22	77.76	97.99	74.22	69.92	88.84	120.00	213.36	48.39
Lipophilicity
iLOGP	3.02	2.91	2.64	2.32	1.75	1.55	3.29	2.41	0.39	4.22	3.40	3.58
XLOGP3	3.47	3.89	3.28	2.37	0.27	-0.96	2.63	4.06	-0.56	5.92	1.91	3.58
WLOGP	2.82	3.26	2.88	2.52	0.96	0.01	3.19	3.74	-0.57	3.57	2.21	3.59
MLOGP	1.55	1.78	1.96	1.06	0.80	0.07	1.31	2.35	-1.30	2.93	0.18	2.35
SILICOS-IT	3.02	3.87	3.23	2.91	1.81	0.93	4.03	3.55	0.69	6.02	-0.55	3.73
Consensus logPo/w	2.78	3.14	2.80	2.24	1.12	0.32	2.89	3.22	-0.27	4.53	1.53	3.77
Water Solubility
ESOL	MS	MS	MS	S	VS	VS	S	MS	VS	PS	MS	S
ALI	MS	MS	MS	S	VS	VS	S	MS	VS	PS	PS	MS
SILICOS-IT	MS	MS	MS	MS	S	S	MS	MS	S	IS	MS	PS
Pharmacokinetics
GI Absorption	High	High	High	High	High	High	High	High	High	High	Low	High
BBB Permeation	Yes	Yes	Yes	No	No	No	Yes	Yes	No	No	No	Yes
Skin permeation, logKp (cm/s)	-5.90	-5.70	-5.80	-6.53	-7.73	-8.72	-6.70	-5.15	-7.66	-5.93	-8.62	-5.81

Fav: Favipiravir, Lop: Lopinavir, Rem: Remdesivir, Hyd: Hydroxychloroquine, TPSA: Topological Polar Surface Area, GI: Gastrointestinal, BBB: Blood Brain Barrier, IS: Insoluble, PS: Poorly soluble, MS: Moderately soluble, S: Soluble, VS: Very soluble.

Table 3

Results of drug-likeness and medicinal chemistry analyses.

Compound	1	2	3	4	5	6	7	8	Fav	Lop	Rem	Hyd
Drug-likeness
Lipinksi	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes	No	Yes
Ghose	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes	No	No	No	Yes
Veber	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes	No	No	Yes
Egan	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes	No	Yes
Muegge	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes	No	No	No	Yes
ABS	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.17	0.55
Medicinal Chemistry
PAINS (alerts)	0	0	0	0	0	0	0	0	0	0	0	0
Brenk (alerts)	0	1	0	0	0	0	0	1	0	0	1	0
Leadlikeness	Yes	No (2)	Yes	Yes	Yes	Yes	No (1)	No (1)	No (1)	No (3)	No (2)	No (2)
SAS	4.27	4.12	3.54	3.53	4.19	4.55	4.03	1.60	2.08	5.67	6.33	2.82

Fav: Favipiravir, Lop: Lopinavir, Rem: Remdesivir, Hyd: Hydroxychloroquine, ABS: Abbott Bioavailability Score, PAINS: Pan Assay Interference Structures, SAS: Synthetic Accessibility Score.

In the study it was predicted that gastrointestinal absorptions of the investigated compounds are high and most of the molecules (except 4, 5 and 6) can cross the blood brain barrier. It was also observed that predicted skin permeations [44] vary between -5.15 and -8.72 cm/s. Results showed that all molecules obey Lipinski, Ghose, Veber, Egan and Muegge rules, and Abbott Bioavailability Scores [45] were predicted to be 0.55 for all molecules.

No alert was observed for any of the investigated compounds in PAINS [46] analyses. In Brenk [47] analyses no alerts was observed for the investigated compounds, except compounds 2 and 8. For compounds 2 and 8, only one alert was observed for each due to bearing isolated alkene group and polycyclic aromatic hydrocarbon group, respectively. In leadlikeness [48] analyses for compound 2, two violations were observed due to its higher molecular weight than 350 g/mol and due to its higher XLOGP3 value than 3.5. For compound 7, one violation was observed due to its higher molecular weight than 350 g/mol and for compound 8, one violation was observed due to its higher XLOGP3 value than 3.5. Finally, results showed that synthetic accessibility scores were obtained in the range of 1.60-4.55.

Conclusion

In this study, eight natural compounds which have recently entered the literature have been investigated computationally for their potential use against SARS-CoV-2 via inhibiting the SARS-CoV-2 main protease. In the study, DFT calculations, molecular docking calculations, molecular dynamics simulations, binding free energy calculations and finally druglikeness analyses and ADME predictions have been performed on the investigated compounds. Results showed that these natural compounds, especially two of them, showed good binding affinity to SARS-CoV-2 main protease and all the molecules obeyed the Lipinksi, Ghose, Veber, Egan and Muegge rules and can be promising structures in the treatment of COVID-19.

Funding

This work was supported by Kocaeli University Scientific Research Projects Unit (Project No: 2011/062). The computational studies reported in this paper were partially performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources) and Kocaeli University.

Availability of data and material

Not applicable

Code availability

Not applicable

Author's statement

Not applicable.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.