Phytochemicals as potential inhibitors for COVID-19 revealed by molecular docking, molecular dynamic simulation and DFT studies.

The COVID-19 pandemic outbreak demands the designing of potential drugs as there is no specific treatment available. Thanks to their safety and effectiveness, phytochemicals have been used to treat various diseases, including antiviral therapeutics. Molecular docking is a simple, quick, and effective way to screen a variety of molecules for structure-based drug design. Here, we investigate molecular docking experiments on compounds present in plant species, Cocculus hirsutus and Rhodiola rosea and show their potential for the treatment of COVID-19. Almost all the components showed higher binding affinity than the built-in ligand, and those with significantly higher binding affinity were explored further. Molecular mechanics-based generalized born surface area calculations were used to re-rank the top candidates, rhodionidin and cocsoline, and their stability in association with viral protease was confirmed. Density functional theory was used for detailed investigations of the geometries, and electrical properties of rhodionidin and cocsoline. Using the frontier molecular orbitals method, the charge transfer within the molecule was calculated. Chemical reactivity and intermolecular interactions were studied using molecular electrostatic potential maps. These in silico discoveries will simulate the identification of powerful COVID-19 inhibitors, and similar research is likely to make a significant contribution to antiviral drug discovery. Supplementary information: The online version contains supplementary material available at 10.1007/s11224-022-01982-4.

Introduction

Infectious diseases pose major threat to public health and social well-being. The coronavirus family gained immense attention, when a new human coronavirus called COVID-19 first appeared in the city of Wuhan, China in December 2019 [1]. The global epidemic of severe acute respiratory syndrome (SARS-CoV-2) has impacted the whole world economy as well[2]. Coronavirus comes under the Coronaviridae [3] family in the order nidovirale [4]. Coronavirus contains crown-like spikes on the virus’ outer surface[5]. Alpha, Beta, Gamma, and Delta coronaviruses are genetically evolved strains of the reference 2019 SARS-CoV-2 genome and are classified as Variants by WHO [6]. The nucleotide sequence present in SARS-CoV-2 indicates that it belongs to the betacoronavirus category [7, 8]. Symptoms of COVID-19 include dry cough, fever, malaise, shortness of breath, and breathing difficullty [9, 10]. Drugs such as hydroxychloroquine [11], choloroquine phosphate [12], lopinavir [13], and remdesivir [14] are the antiviral therapies available for the coronavirus as of now. However, their effectiveness remains uncertain and needs further evaluation. Several therapeutic advancements are underway in order to find an effective medication that can directly target viral proteins.

Plant species are immense source of biologically important secondary metabolites such as flavonoids, alkaloids, glycosides, tannins, terpenoids, polyphenols, and saponins [15, 16]. Considering the multiple biological activity, researchers are focusing on the several plant-based compounds that can act as potent drug components for various diseases including SARS-CoV-2. However, detection of potential inhibitors using traditional methods is expensive as well as time-consuming. Therefore, in recent years, the use of in silico identification of druggable molecules is gaining tremendous attention as they offer fast and easy detection of potent drugs against various diseases [17, 18]. The discovery and isolation of medicinally important therapeutic compounds from higher plants for various diseases have been well established [19-23].

The screening of herbs that may contain anti-coronavirus compounds from a huge library of species is a big challenge. Cocculus hirsutus L.(C.hirsutus), a climbing shrub abundant with alkaloids, flavanoids, and phenolic compounds, comes under Menispermaceae family [24]. The plant has historically been disparaged for its peculiar property of healing many kinds of cuts and wounds and is common in tropical and subtropical environment [25]. It has been shown that this plant is one of the best candidates for treating various diseases such as gonorrhoea, spermatorrhoea, urinary disorders, diarrhoea, and hyperglycemia [26]. The phytochemicals present in this plant are involved in antimalarial [27], antibacterial [26], hyperglycaemic [28], anti-inflammatory, analgesic [28], diuretic, and laxative [29] activities. The genus Rhodiola, belonging to the family Crassulaceae, is composed of almost 200 species. The well-known species among them is Rhodiola rosea L(R. rosea). This plant is also called golden root or arctic root [30], is mainstream among herbal medicines, and is commonly seen at higher altitudes in the arctic and in mountain ranges throughout Europe and Asia [31]. The medicinal application of R.rosea includes its ability to control psychological stress and mental strength, change the neurotransmitter levels and Central Nervous System (CNS) activity, resistance to high altitude sickness, treat fatigue [32-36], and also act as an anti-depressant and anti-inflammatory drug [37]. Phytochemical analysis on R.rosea showed that the plant is well-equipped with bioactive components like organic acids, flavonoids, monoterpenes, triterpene, and tannins, large amounts of phenolic compounds, and some specific phenylpropane derivatives like rosavins [38].

The key druggable target of SARS-CoV-2 involves 3-chymotrypsin-like protease (3CLpro) [39, 40], papain-like protease (PLpro) [41-44], RNA-dependent RNA polymerase [45], and spike (S) proteins [46, 47]. In this work, 6LU7, main protease (Mpro) of COVID-19, crystallized by Liu et al. (2020), is used as a potential protein target for the phytochemicals selected. The in silico docking model is the main protease, which regulates the major functions of the virus and has a strongly preserved catalytic domain from the SARS virus [48]. Some of its roles involve the virus’ replication, making it an attractive target for drug development [49]. Usable bioinformatics techniques such as molecular docking can be implemented using comprehensive 3D structures of related biomolecule and phytochemicals to rapidly classify promising and potential drug candidates. In order to completely unravel the fundamental interactions responsible for drug-receptor binding and subsequent structural modifications, this can be accompanied with atomistic molecular dynamics, leading to changes in biological system function. The current work focuses on the identification of bioactive components present in the herbs C.hirsutus and R.rosea and their antiviral activity towards COVID 19 main protease 6LU7 with the aid of molecular docking and molecular dynamic simulation studies. Towards the end we use density functional theory (DFT) at the B3LYP/6–31 +  + G(d,p) level of theory to find the optimal structures of the molecules from each plant that have the maximum binding affinity with the specified target. Following that, their reactivities were predicted at the same level of theory using frontier orbital investigations. Furthermore, molecular electrostatic potential surfaces were investigated to determine which parts of a molecule are the most reactive nucleophilic and electrophilic against reactive biological potentials.

Materials and methods

Molecular docking

Preparation of protein/receptor

The crystal structure of COVID-19 main protease in complex with an inhibitor N3 (PDB ID- 6LU7) with resolution 2.16 Å was downloaded from RCSB protein data bank database (www.rcsb.org) as.pdb format. The 6LU7 protein contains two chains, A and C, in which A represents the protein and C corresponds to the inbuilt ligand. Prior to the detailed docking analysis, chain A was chosen as the target for this study and then the preparation of the protein was performed using AutoDock tool. The selected protein complex contains an inbuilt ligand 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.

Both the water molecules and non-interacting ions, along with the native ligand, were eliminated from the crystal structure. The missing hydrogens were added in order to alleviate the stress of the crystal structure and make the protein accessible for use in the AutoDock docking simulation program. After the structural minimization, the protein was prepared using AutoDock Tools (version1.5.6) (ADT), graphical user interface, which involves the addition of hydrogen atoms, Gasteiger charge calculation, and merging of non-polar hydrogens to carbon atoms. The generated macromolecular structure was then saved as pdbqt file.

Ligand preparation

Most of the ligands were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The downloaded ligands in.sdf format were converted into.pdb file using the Open Babel software (version 2.4.1). The ligand structures that were unavailable in the PubChem database were created using Maestro and the structures were downloaded in pdb format. Ligand input files for docking was prepared using AutoDock Tools and saved as pdbqt files.

Docking

For the detailed analysis of intermolecular interactions and binding modes between the protein and specific phytochemicals present in C. hirsutus and R. rosea, the molecular docking studies were executed for selected ligands using the AutoDock Vina software. A grid box with dimension X:15, Y:22, Z:20 Å with a grid spacing of 1.0 Å centered at X: −9.491, Y:11.784, Z: 65.363A0 was identified as protein target docking site. The interactions between active sites in the target and ligand conformation, including the type of interaction and bond distances, were identified using Discovery Studio Visualizer.

Binding energy calculations using prime/MMGBSA studies

The calculation of free energy provides a quantitative estimate of protein–ligand interactions that can predict the stability of protein–ligand complexes. Molecular mechanics-generalized born surface area (MM-GB/SA) was conducted with Prime module of Schrodinger 2018 to calculate binding energy of protein–ligand complexes. MM-GBSA method in Prime was used for rescoring the docked pose of ligands together with calculation of coulomb energy, covalent binding energy, hydrogen-bonding correction, pi-pi packing correction, lipophilic energy, generalized born electrostatic solvation energy, and Van der Waals energy.

Molecular dynamics simulation

To forecast the stability of protein–ligand complexes, molecular dynamic simulations were performed in Desmond package of the Schrodinger suit. The chosen protein–ligand complexes were first immersed in the TIP4P water box, extending 10 Å beyond any atoms of the complex. The orthorhombic boundary condition box is selected, and buffer method is used for box size calculation. To neutralize charges, counter ions (3 Na+) were added. The MD was conducted at a temperature of 310 K and 1.63 bar pressure over 100 ns in the NPT ensemble with recording intervals of 12 ps for trajectory. Simulations were performed with the force field of OPLS-3e. The Desmond simulation interaction diagram method from Maestro was used to sketch plots and figures.

Quantum chemical analysis

Density functional theory (DFT) is an intriguing method for probing many biological features of interest. The current study employs DFT analysis to highlight the electronic features of the phytochemicals. The geometry optimizations of the selected compounds cocsoline (from C.hirsutus) and rhodionidin (from R.rosea) were carried out using Gaussian16 at the B3LYP/6–31 +  + G(d,p) level. Based on the optimized structures in the gas phase, electronic characteristics such as HOMO–LUMO energies and molecular electrostatic potential (MESP) were computed. The natural population analysis, as well as global reactivity descriptors, was further investigated for the selected compounds.

Results and discussion

This section is organized as follows. First, we identified the chemical components of C.hirsutus and R.rosea and the identified components were docked with 6LU7, the SARS-CoV-2 target protein. The components with highest binding affinity obtained using docking were re-ranked using MMGBSA calculation. Then, we evaluated the stability of the complex formed between the components having the highest binding affinity with 6LU7, using RMSD values obtained from the MD simulation analysis. The position of binding modes of components, rhodionidin and cocsoline (with the highest binding affinity with 6LU7) was verified using the MD modelling counterpart. Finally, the density functional theory calculations were performed to investigate the structure and reactivity of top ranked candidates.

Molecular docking analysis of phytochemicals present in C.hirsutus and R.rosea

The essential parts of the structure-based drug designing phase involve the docking of small molecules into the receptor’s active site and evaluating the binding affinity of ligand towards target protein. In the current study, the selected protein, 6LU7 contains a peptide like inbuilt ligand N3 (Fig. 1) that exhibits very potent inhibition against COVID-19 virus Mpro [50]. In order to compare the binding affinity as well as the inhibition potency of selected phytochemicals, the docking process was also performed on the native ligand.

Fig. 1

Two-dimensional image of inbuilt ligand N3, figure taken from reference [51]

The plant C.hirsutus as a rich source of potential anti-coronavirus drugs could be evaluated by analyzing molecular docking results and comparing these with the values of in-built ligand. The structures of the phytochemicals with highest binding affinity values are presented in Fig. 2. The binding affinity values of the phytochemicals (along with PDB ID) present in C.hirsutus as calculated from the software AutoDock Vina are shown in Table S1 and the binding affinity values of the components with high inhibition efficiency than that of in-built ligand N3 are shown in Fig. 3a. The binding affinity of N3 with the target protein obtained from the docking study was −7.0 kcal/mol. The lowest binding energy value of −9.1 kcal/mol was observed for the phytoconstituent, cocsoline. The flavonoids rutin, liquiritin, and quercetin with binding energy values −8.8, −7.5, and −7.5 kcal/mol respectively prove that they could also act as possible inhibitors against 6LU7 protein target. Isotrilobine, coclaurine, cohirsine, and lupeol, exhibiting binding affinity values like −8.6, −7.6, −7.4, −7.3, and −7.1 kcal/mol respectively, are still higher than that of binding affinity of native ligand N3.

Fig. 2

Selection of structures of phytochemicals present in R. rosea that show strong binding affinity (less than or equal to −7.5 kcal/mol) towards SARS-CoV-2

Fig. 3

a The binding affinity values of the components in C.hirsutus with high inhibition efficiency than that of in-built ligand N3 as calculated from the software AutoDock Vina. b The binding affinity values of the components in R.rosea with high inhibition efficiency than that of in-built ligand N3 as calculated from the software AutoDock Vina. *Native ligand N3 shown in red dot

Docking analysis of phytochemicals of the plant C.hirsutus shows that all of them show strong interactions with the surrounding amino acid residues. The two- and three-dimensional interactions of components of C.hirsutus having low binding energy values <  −7.4 kcal/mol as well as the inbuilt ligand N3 with 6LU7 are represented in Table 1. Cocsoline, having the highest binding affinity, possess hydrogen bond interaction between its hydroxyl group and amino acid residues LEU141 and GLY143 near to it. This molecule, within the protein active site, shows hydrophobic interactions such as alkyl and pi-alkyl interactions with CYS145, MET165, and PRO168. The molecule has interacted well with a number of amino acid residues like ASN142, HIS41, THR190, and GLN189 near to it. The structure, 3D images of the phytochemicals within the active site of 6LU7, and 2D interaction of the phytochemicals present in the plant C.hirsutus, with binding energy values greater than −7.4 kcal/mol with the target, are given in Table S3.

Table 1

The 3D images of inbuilt ligand N3, and the phytochemicals of the plant, C.hirsutus having low binding energy values <  −7.4 kcal/mol within the active site of 6LU7 and their interactions with the target protein

The majority of the chemicals found in R.rosea showed good binding affinity values with the specified protein target, 6LU7, according to data from molecular docking. The most of the binding affinity values were found to be greater, with some even being comparable to the native ligand N3. The molecular docking results of bioactive components (along with PDB ID) present in R.rosea with SARS-CoV-2 are shown in the Table S2 and the binding affinity values of the components with high inhibition efficiency than that of in-built ligand N3 are shown in Fig. 3b. The binding affinity of phytochemicals present in R.rosea ranges from −9.6 to −4.0 kcal/mol and most of the components show higher inhibition tendency towards SARS-CoV-2 compared to the native ligand. Among the bioactive components selected, rhodionidin, rhodiolgidin, rhodalin, rhodalidin, rhodionin, and rhodiolgin, coming under the flavonoid category exhibited lowest binding energy values such as −9.8, −8.9, −8.6, −8.5, −8.0, and −8.0 kcal/mol respectively. The structures of phytochemicals of R.rosea having low binding energy values less than or equal to −8.0 kcal/mol are shown in Fig. 4. The flavonoids such as kaempferol, rhodiosin, gossypetin, and herbacetin also show high binding affinity values which are still higher than the binding affinity of inbuilt ligand (−7.0 kcal/mol). Rosavin, rosarin, and sachaliside, belonging to the phenyl proponoid group, also showed a high affinity towards the target protein. Docking studies revealed that rosin and vimalin (belonging to phenylproponoids) and rhodioloside (belonging to phenylethanoids) show a binding energy value with the SARS-CoV-2 which is equal to the binding energy of inbuilt ligand with the same target. Other components like picein, rosiridin, gallic acid, rosiridol, geraniol, cinnamyl alcohol, limonene, alpha-pinene, and decanol exhibit low inhibition for the target. More interestingly, these values are still comparable with the native ligand.

Fig. 4

Selection of structures of phytochemicals present in R. rosea that show strong binding affinity (binding energy values less than or equal to -8 kcal/mol) towards the SARS-CoV-2

The three- and two-dimensional binding interactions of best scoring components present in R.rosea having low binding energy values <  −7.9 kcal/mol with 6LU7 are provided in Table 2. The major interactions involved are conventional hydrogen bonding, alky, pi- alkyl, pi-pi, pi-S, pi-cation, carbon-hydrogen bond, and pi donor hydrogen bond and all of the chemical components shows very strong interaction with the protein studied.

Table 2

The 3D images of the phytochemicals present in the plant R. Rosea having low binding energy values <  −7.9 kcal/mol within the active site of 6LU7 and their interactions with 6LU7

The detailed visualization of stable conformer of rhodionidin, the top-ranked (binding energy value of −9.6 kcal/mol) SARS-CoV-2 inhibitor at the binding site of 6LU7, reveals the existence of a strong hydrogen bond interaction between hydroxyl groups present in the ligand and GLY143, SER144, LEU141, and HIS163 amino acid residues of the target protein. A π-π interaction exists between the inhibitor and HIS41 residue present in the active site of 6LU7. Alkyl and aromatic fragments present in the inhibitor form hydrophobic interactions such as alkyl and pi alkyl interaction with LEU27, MET49, PRO168, and MET165 residues. Furthermore, the molecule is well fixed within the cavity of the active site.

From the detailed interaction studies of phytochemicals present in the herb R. rosea showing high binding energy with selected protein proved that all components are potential inhibitors and possess strong interactions with the target. The structure, 3D images of the phytochemicals within the active site of 6LU7, and 2D interaction of the phytochemicals present in R. rosea, with binding energy values less than -8.0kcal/mol with the target, are given in Table S4.

Molecular dynamic simulation using Desmond package

Since docking is a static view of the molecule’s binding position in the protein’s active site, molecular dynamics (MD) simulation attempts to quantify atomic motions over time by integrating the classical motion equation of Newton. The dynamic behavior of the molecular system is simulated in MD, in order to assess the stability of the protein–ligand complex. For the MD analysis with the OPLS3e force field, the docked conformers of cocsoline and rhodionidin with the highest predictive binding energy of −9.1 and −9.6 kcal/mol were therefore used. The dynamic features of these two selected systems were calculated during the period of 100 ns simulation. We evaluated the stability of the docked compounds in the active pocket of 6LU7 by using root mean square deviation (RMSD) and its effects on the stability of the entire system.

RMSD plot in Fig. 5a indicates the rhodionidin:6lu7 complex MD trajectory for 100 ns. The complex appears to be stabilized within a short time with respect to the reference frame at time 0 ns, during the simulation phase. The RMSD variation of 6lu7 was only around 1.6 Å which can be considered insignificant. This in turn means that the protein RMSD appears to be stable during the simulation, indicating that the protein has not undergone major conformation changes. Fluctuation within 20–40 ns of the trajectory was observed in the case of ligand RMSD which indicate major conformational changes of the ligand. After 40 ns, the ligand RMSD stabilized again and remained constant until the end of the simulation period of 100 ns. The ligand has undergone conformational changes with an overall RMSD of 6 Å, possibly because of a large number of rotatable bonds that contribute to its flexibility. The protein RMSD appears to be constant with negligible variation which means that the protein is undergoing minimum conformational change during the simulation. The rhodionidin:6lu7 complex’s overall RMSD plot shows that the ligand is less stably connected to the protease binding site.

Fig. 5

RMSD analysis for MD simulation trajectories for a rhodionidin:6lu7 complex and b cocsoline:6lu7 complex

The RMSD plot in Fig. 5b shows that the cocsoline:6lu7 complex stabilized immediately after the simulation started, relative to the reference frame at 0 ns. During the simulation, the complex appears to be stabilized with regard to the reference frame at time 0 ns. It is possible to see a small divergence at the end of the simulation, i.e., 80 ns. The fluctuation, however, lies within the acceptable range of 1–3 Å, and can therefore be considered non-significant. In this complex, since the ligand and protein backbone RMSD plots were lying over each other, it is possible to conclude a stable complex formation. RMSD value analysis reveals that the complex cocsoline:6lu7 was more stable compared to the rhodionidin:6lu7 complex.

The position of the binding modes of compounds rhodionidin and cocsoline has also been verified by the MD modelling counterpart. During the simulation, protein interactions with the ligand were regulated. Such interactions can be arranged and analyzed as depicted in Fig. 6. Protein–ligand interactions have been classified as hydrogen bonds, hydrophobic, ionic, and water bridges. In order to keep the protein–ligand complex stable, the form of amino acid residues present in the active site of the target protein plays a crucial role. Figure 6a indicates the type of contacts maintained in the entire trajectory of rhodionidin:6LU7 complex. In the case of rhodionidin, GLN189, GLU166, and GLN192 were able to hold down the hydrogen bonding interactions during 85%, 48%, and 47% of the time, respectively. It also showed a greater number of hydrogen bonding contacts with THR26, SER46, HIS164, and THR190 around 10–20% of the time. Hydrophobic interactions were found with HIS41, MET165, CYS145, and MET49, during 60%, 55%, 15%, and 10% respectively. Water bridge hydrogen bonding, where one water molecule is used as a bridge between ligand and protein, was also formed with GLU166, GLN189, THR24, THR25, THR26, SER46, ASN142, ARG188, THR190, and GLN192.

Fig. 6

a Rhodionidin: 6LU7, b cocsoline: 6LU7 contact throughout the trajectory

The various types of interaction between cocsoline and 6LU7 at the active site during the entire simulation period are represented in Fig. 6b. Cocsoline showed hydrogen bonding interaction with THR26 and HIS41 which were 49% and 30% of the time, respectively. It also exhibits weak hydrogen bonding contacts with GLN189, GLY143, HIS164, and GLU166 less than 10% of the time. Cocsoline possesses multiple water bridge hydrogen bonding with THR190, GLN192, GLN189, GLU166, CYS145, GLY143, SER144, TYR54, ASN119, THR25, HIS41, MET49, ASN142, HIS163, HIS164, ASP189, and ARG188 residues.

Relevant timeline contacts of the compounds rhodionidin and cocsoline with the amino acid residues that were present in the targets are also studied and provided in Fig. 7. The darker color reflects the high number of ties to the amino acid.

Fig. 7

a Rhodionidin: 6LU7 and b cocsoline: 6LU7 interactions shown in each trajectory frame by the active site amino acids; no interactions are represented by white while stronger interactions by the dark color

MMGBSA analysis

MMGBSA is a popular analysis method to calculate the average binding energy of equilibrated MD trajectory of the systems. MMGBSA calculations were conducted in order to provide a better ranking of the ligands, and for the determination of predictive binding energies. The binding free energies of the two selected systems indicate the binding affinity of ligands with 6LU7. The values obtained after measurement were estimated free binding energies with a more negative value, suggesting a stronger binding energy. Table 3 depicts the average binding energy, coulomb energy, covalent binding energy, hydrogen-bonding correction, pi-pi packing correction, lipophilic energy, generalized born electrostatic solvation energy, Van der Waals energy of rhodionidin, and cocsoline with 6LU7. Rhodionidin exhibited the highest negative binding energy of −80.40 kcal/mol. MMGBSA calculation revealed that cocsoline showed an approximate binding energy value of −75.07 kcal/mol.

Table 3

MMGBSA analysis of rhodionidin and cocsoline with 6lu7

Molecule	∆G Binding(kcal/mol)	Coulomb(kcal/mol)	Covalent(kcal/mol)	H-bonding(kcal/mol)	Bind packing(kcal/mol)	Lipo(kcal/mol)	Solv_GB(kcal/mol)	vdW(kcal/mol)
Rhodionidin	−80.40	−41.28	7.53	−3.91	−3.05	−16.78	39.04	−61.95
Cocsoline	−75.07	−18.38	3.82	−1.49	−1.92	−24.17	29.18	−62.10

DFT studies

Geometry optimization

The optimized structures of cocsoline and rhodionidin obtained at B3LYP/6–31 +  + G(d,p) level are given in Fig. 8, and the corresponding geometrical parameters are given in Tables S5–S10. The total energy of cocsoline and rhodionidin computed at the same basis set are −1800.92 and −2251.14 Hartree respectively. A molecule’s dipole moment is a three-dimensional vector that depicts the molecular charge distribution. As a result, it can be used as a descriptor to characterize the charge flow within a molecule. The DFT/ B3LYP/6–31 +  + G(d,p) computations showed that, the rhodionidin and cocsoline possessed a dipole moment value of 8.88 Debye, and 2.54 Debye respectively. Rhodionidin has a greater dipole moment than cocsoline, indicating that it is softer. The large number of hydroxyl groups in the structure of rhodionidin may facilitate a greater number of hydrogen bond formation.

Fig. 8

Optimized geometry of a cocsoline and b rhodionidin obtained from DFT/ B3LYP/6–31 +  + G(d,p) calculation

Natural population analysis and molecular electrostatic potential calculations

Molecular polarizability, dipole moment, electronic structure, molecular reactivity, and many other aspects of molecular systems are influenced by atomic charges. The creation of donor and acceptor pairs involving the charge transferring molecule is suggested by the charge distributions over the atoms. The distribution of electrons in various sub shells of atomic orbitals is revealed by a natural population analysis (NPA) of an organic molecule [52]. The natural charges on atoms of cocsoline and rhodionidin are given in Tables 4 and 5 respectively. The value of natural charges informs that the atom O10 and O5 were showed more electronegativity in cocsoline and rhodionidin respectively than all other oxygen atoms. Among the two nitrogen atoms present in rhodionidin, N7 was more electronegative than N6 atom. On other hand H67 of cocsoline and H73 of rhodionidin were possessed high electropositive value.

Table 4

The natural charges on atoms of cocsoline calculated under DFT/ B3LYP/6–31 +  + G(d,p)

Atom with number	Natural charge	Atom with number	Natural charge	Atom with number	Natural charge
O1	−0.61608	C25	0.11209	H49	0.17917
O2	−0.60443	C26	0.42277	H50	0.16055
O3	−0.56579	C27	−0.04182	H51	0.20232
O4	−0.56832	C28	−0.58562	H52	0.18214
O5	−0.74962	C29	0.31144	H53	0.18942
O6	−0.74207	C30	0.30742	H54	0.18801
O7	−0.73267	C31	0.31173	H55	0.19396
O8	−0.75343	C32	−0.33928	H56	0.22399
O9	−0.74592	C33	−0.26458	H57	0.20485
O10	−0.77376	C34	0.40427	H58	0.20846
O11	−0.75721	C35	0.45172	H59	0.46319
O12	−0.48107	C36	0.32245	H60	0.47688
O13	−0.66685	C37	0.19105	H61	0.47221
O14	−0.68832	C38	−0.13628	H62	0.47816
O15	−0.68158	C39	−0.13762	H63	0.47442
O16	−0.66728	C40	−0.14611	H64	0.49738
C17	0.07465	C41	−0.26048	H65	0.24829
C18	0.08538	C42	−0.28936	H66	0.4742
C19	0.07487	C43	0.33181	H67	0.50614
C20	0.09108	H44	0.20341	H68	0.23775
C21	0.08571	H45	0.16869	H69	0.24768
C22	0.0825	H46	0.179	H70	0.50011
C23	0.43941	H47	0.18083	H71	0.22267
C24	0.06096	H48	0.20121	H72	0.20166
				H73	0.47060

Table 5

The natural charges on atoms of rhodionidin calculated under DFT/ B3LYP/6 −31 +  + G(d,p)

Atom with number	Natural charge	Atom with number	Natural charge	Atom with number	Natural charge
O1	−0.52301	C 25	0.25684	H49	0.18135
O2	−0.51843	C 26	−0.12553	H50	0.20602
O3	−0.59358	C 27	0.43116	H51	0.21479
O4	−0.52784	C28	−0.0429	H52	0.21451
O5	−0.66359	C29	−0.21052	H53	0.16073
N6	−0.58127	C30	−0.04229	H54	0.19579
N7	−0.67442	C31	−0.18599	H55	0.35232
C8	0.10144	C32	−0.1904	H56	0.19146
C9	−0.12319	C33	−0.25236	H57	0.19031
C10	−0.0301	C34	−0.21128	H58	0.16462
C11	−0.18439	C35	−0.23911	H59	0.21806
C12	−0.38164	C36	−0.23511	H60	0.21934
C13	−0.41174	C37	0.28812	H61	0.21029
C14	0.27304	C38	0.2822	H62	0.2174
C15	−0.0229	C39	−0.26334	H63	0.2063
C16	−0.0585	C40	0.27411	H64	0.21103
C17	−0.03401	C41	−0.20129	H65	0.2162
C18	0.1735	H42	0.21273	H66	0.20104
C19	−0.40005	H43	0.19657	H67	0.22018
C20	−0.27856	H44	0.18965	H68	0.22056
C21	−0.17441	H45	0.20447	H69	0.20191
C22	−0.35356	H46	0.20965	H70	0.18408
C23	0.22297	H47	0.20441	H71	0.16365
C24	−0.40195	H48	0.22366	H72	0.17909
				H73	0.46795

The three-dimensional charge distributions of molecules will be depicted by the molecular electrostatic potential surface. The molecular electrostatic potential (MEP) surface can also be viewed as a map of electron excess and deficient areas [53, 54]. The MEP must be determined to validate the evidence of the drug’s reactivity as inhibitors, despite the fact that the MEP indicates the molecule size and shape of the positive, negative, and neutral electrostatic potentials. These could be used to anticipate physicochemical property relationships based on the molecular structure of medications in development. Furthermore, the molecular electrostatic potential can be used to evaluate a drug’s responsiveness to electrophilic and nucleophilic assaults. The MEP diagram of both cocsoline and rhodionidin calculated at DFT/ B3LYP/6–31 +  + G(d,p) is given in Fig. 9. In the MEP, the largest negative region, which is shown red, is the favored site for electrophilic attack. As a result, an attacking electrophile will be drawn to the negatively charged sites, whereas the blue regions will attract the reverse. The change in the drug’s binding affinity with the active site receptor could be due to differences in the mapping of the electrostatic potential around it.

Fig. 9

MEP of a cocsoline and b rhodionidin under DFT/B3LYP/6–31 +  + G(d,p)

Frontier molecular orbital analysis

The value of highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) energy can be used to determine a molecule’s ability to donate and receive electrons. The fundamental idea of frontier molecular orbital (FMO) theory can be summarized in the form of a simple rule that states that the maximum positive overlap between LUMO (empty state) and HOMO (filled state) orbitals is required for a simple reaction course. HOMO (highest occupied molecular orbital) is related to ionization potential, while LUMO (lowest unoccupied molecular orbital) is related to electron affinity [55]. These molecular orbitals are important in electrical and optical characteristics, and pharmacological investigations, as well as providing biological mechanism information [56]. The energy gap of the frontier molecular orbitals (i.e., FMO) supports the stability of structure. FMOs also provide information on a molecule’s kinetic stability and chemical reactivity. The FMOs theory demonstrated that the HOMO and LUMO energy levels had the greatest impact on the bioactivities of small structural medicines. The HOMOs are the ones which provide electrons, whereas the LUMOs are the ones which accept them. Furthermore, the FMOs aid in the prediction of a molecule’s most reactive site [57]. The FMOs in the electronic transitions and their energies difference (Fig. 10) are determined in order to anticipate the energetic behaviors and reactivity of chloroquine and chloroquine phosphate against COVID-19 virus. In Fig. 10, the negative phase is represented by the color green, whereas the positive phase is represented by the color red.

Fig. 10

Frontier molecular orbitals of a cocsoline and b rhodionidin obtained from DFT/ B3LYP/6–31 +  + G(d,p) calculations

Rhodionidin showed the most lying HOMO than cocsoline, and consequently, it could be a better electron donor drug. The energy gap of rhodionidin and cocsoline is 3.6662 and 4.7775 eV respectively. Rhodionidin’s small energy gap allows electrons to move more freely, making the molecule soft and reactive and this result is in strong agreement with the high dipole moment of rhodionidin. The HOMO of a certain drug and the LUMO with the adjacent residues could share the orbital interactions during the binding process.

Chemical reactivity descriptors

The EHOMO and ELUMO are markers that can be used to forecast a molecule’s ionization potential (I =  −EHOMO) and electron affinity (A =  −ELUMO). Other chemical reactivity descriptors such as hardness (η), softness (S), electronegativity (χ), chemical potential (μ), and electrophilicity index (ω) are estimated using the frontier molecular orbitals and are shown in Fig. 11.

Fig. 11

Global descriptive parameters of a cocsoline and b rhodionidin as obtained from DFT/ B3LYP/6–31 +  + G(d,p) calculations

The electronegativity (χ = (I + A)1/2) value predicts the molecule’s ability to attract electrons, i.e., Lewis acid, whereas lower values of χ indicate a suitable base. The χ value is found to be 3.15 and 4.09 for cocsoline and rhodionidin, respectively. The global hardness (η = (I −A)1/2) is a measure of their charge transfer prohibition, while the global softness (S = 1/2η) characterizes a molecule’s ability to accept electrons. Because they may easily transmit electrons to acceptors, soft molecules have a smaller energy gap between frontier molecular orbitals and are more reactive than harder molecules. The hardness of cocsoline (2.39 eV) was higher than that of rhodionidin (1.83 eV) whereas the softness of cocsoline (0.21 eV) is lower than that of rhodionidin (0.27 eV). The electrophilicity index (ω) = μ2/2η determined for cocsoline and rhodionidine is 2.07 and 4.56 eV, respectively. The chemical potential (μ =  −χ) is found to be −3.15 eV for cocsoline and −4.09 eV for rhodionidin.

Conclusion

Components present in the two medicinal herbs R.rosea and C.hirsutus were subjected to molecular docking studies in order to identify possible COVID-19 Mpro inhibitors. Using the AutoDock Vina software, the binding energy values of all the components present in both plants were determined. The findings of in silico molecular docking studies of R.rosea and C.hirsutus indicate that most of the phytochemicals have moderate to potent antiviral activity against SARS-CoV-2 Mpro. The components’ high affinity is also found to be dependent on the nature and strength of bonding within the protein’s active site. The flavonoid rhodionidin, which is found in the R.rosea plant and has the lowest binding energy of −9.6 kcal/mol, and the alkaloid cocsoline, which is found in the species C.hirsutus and has a binding energy of −9.1 kcal/mol, have been described as the best potential drug against the novel coronaviruses. MMGBSA analysis was used to help re-rank these ligands showing the highest binding affinity for 6LU7. MD simulations of rhodionidin and cocsoline with the target indicate that the complex cocsoline:6lu7 was more stable compared to the rhodionidin:6lu7 complex. The DFT/ B3LYP/6–31 +  + G(d,p) approach was used to optimize the molecular structures of cocsoline and rhodionidin (having highest binding affinity with the target protein), as well as to identify their geometrical characteristics. Frontier orbitals, gap energies, natural population analysis, and reactivity descriptors are some of the molecular features that have been discussed. According to the findings, the minimal energy gap is associated with a high binding affinity for rhodionidin. The molecule gets softer and more reactive as the gap energy decreases, making electron flow easier. Following that, the estimated MEP maps reveal that positive potential sites are more advantageous for nucleophilic attack, while negative potential sites are more favorable for electrophilic attack. These findings further motivate in vitro and in vivo studies of the identified components for the treatment of COVID-19. These findings are also expected to stimulate further research into the production of safe and successful anti-coronavirus or other antiviral drugs derived from naturally occurring compounds.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOC 4726 KB)