Computational studies of Bridelia retusa phytochemicals for the identification of promising molecules with inhibitory potential against the spike protein and papain-like protease of SARS-CoV-2

1.            Introduction

Severe Acute Respiratory Syndrome (SARS), a viral respiratory illness, is caused by a coronavirus called SARS-associated coronavirus (SARS-CoV). SARS-CoV-2 is the virus responsible for the dreadful infectious coronavirus disease 2019 (COVID-19). The symptoms of COVID-19 are fever, cough, tiredness, loss of taste and smell, headache, sore throat, aches & pains, diarrhea, difficulty in breathing, and chest pain. Majority of the infected people who developed mild to moderate symptoms can recover without hospitalization whereas medical attention is required in case of serious symptoms [1]. By the end of 2019, a novel coronavirus designated as SARS-CoV-2 emerged in the city of Wuhan, China, and made an outbreak of unusual viral pneumonia [2].

To date, a total of 539,119,771 confirmed cases of COVID-19, of which 6,322,311 deaths have been reported by the World Health Organisation (WHO) [3]. The SARS-CoV-2 is a member of the beta coronavirus (beta-CoV) family, which includes RNA viruses with spikes resembling crowns on the surface of the coronavirus particles, similar to the SARS (severe acute respiratory syndrome) and MERS (Middle East Respiratory Syndrome) viruses. However, the WHO estimates that the death rate for SARS-CoV-2 is between 2 and 3 percent, compared to 14 and 35 percent for SARS and MERS, respectively. However, the SARS-CoV-2 possesses characteristics such as rapid person-to-person transmission, asymptomatic transmission, protracted symptom development, as well as significantly higher fatality rates in the older population [4].

The pandemic has increased the demand for critical care, placing enormous strain on the healthcare systems of many nations. Many therapies were suggested for the treatment of SARS-CoV-2 during the peak COVID-19 pandemic. The WHO has launched a trial, SOLIDARITY, to focus on testing the four most promising COVID-19 treatments namely, remdesivir, chloroquine and hydroxy-chloroquine, lopinavir, and lopinavir and ritonavir and interferon-beta. Chloroquine/hydroxy-chloroquine and lopinavir/ritonavir were removed from the COVID-19 treatment protocols in June 2020 due to possible risks and uncertainty of their benefits. Remdesivir was the first and only drug approved by the FDA on 22 October 2020 [5,6].

However, SARS-CoV-2 kept mutating leading to the emergence of deadlier variants with higher transmissibility which might weaken the efficacy of the available vaccines and antiviral medications. Therefore, the search for novel or adjuvant anti-SARS-CoV-2 medications is urgently required. Investigating the bioactive compounds of plants appears to be a promising approach for the identification and development of new medications for COVID-19 [7]. Numerous plants have been reported to possess antiviral phytoconstituents that may be beneficial against SARS-CoV-2 namely, Lycoris radiata, Artemisia annua, Pyrrosia lingaua, Lindera aggregate, Isatis indigotica, Torreya nucifera, Houttuynia cordata, Curcuma longa, and Curcuma xanthorrhiza [8,9].

Bridelia retusa (L.) A. Juss. is a deciduous tree found throughout India that grows to be about 20 meters tall. The plant is used as a traditional herbal remedy. Natives use the stem bark and roots to treat rheumatism and as astringent agents [10]. During fever, a powdered stem bark mixed with water is administered [11]. The plant is said to have a variety of biological activities such as anti-viral, respiratory, cardiovascular, and anticancer activities [12]. In-vivo studies have also revealed that B. retusa has antinociceptive, anti-inflammatory, antimicrobial, and anti-fungal properties and also stimulates cell-mediated immunity [13-15]. The facts presented above supports the rationale for investigating B. retusa phytocompounds as a potential antiviral agent against SARS-CoV-2. Phytoconstituents like tannins, triterpene ketone, decanoic acid octadecyl ester, stigmasterol, and dehydrostigmasterol are found in the bark. The fruits contain β-sitosterol, gallic acid, and ellagic acid and the leaves contain crude proteins [16].

In our present study, we have selected 14 phytocompounds namely lupeol, friedelin, gallic acid, ellagic acid,  (E)-4-(1,5-dimethyl-3-oxo-1,4-hexadienyl) benzoic acid, (E)-4-(1,5-dimethyl-3-oxo-1-hexenyl) benzoic acid, (R)-4-(1,5-dimethyl-3-oxo-4-hexenyl) benzoic acid, (-)-isochaminic acid, (R)-4-(1,5-dimethyl-3-oxohexyl) benzoic acid (ar-todomatuic acid), 5-allyl-1,2,3-trimethoxybenzene (elemicin), (+)-sesamin, β-sitosterol, stigmasterol, and 4-isopropylbenzoic acid (cumic acid) that have been reported to be present in B. retusa and were investigated for their inhibitory potential against the main protease of SARS-CoV-2 through in-silico approach [17]. In SARS-CoV-2, the spike protein (S) help the virus to penetrate the host cell [18] and the papain-like protease (PL pro) protein is necessary for the processing of viral polyproteins to generate a functional replicase complex and enable the spread of the virus [19]. In the present study, the inhibitory potential of the selected 14 phytocompounds against these protein targets will be investigated with in-silico techniques.

2.            Materials and Methods

2.1.        Retrieval of Target protein

The X-ray crystal structure of SARS-CoV-2 S-protein (PDB ID: 6M0J) and PL pro (PDB ID: 7OFT) are selected as the target proteins. They were downloaded in ‘.pdb’ format from the RCSB-PDB website. The crystal structures of the target chains of the proteins are shown in Figure 1 and Figure 2. The S-protein was found to have two Chains (Chain A and Chain E) and the PL pro has only one Chain (Chain A). The co-crystal inhibitor of the S-protein (2-acetamido-2-deoxy-beta-D-glucopyranose, NAG) and PL pro protein (p-hydroxybenzaldehyde, HBA) was identified and they manually prepared with the Chem Draw Professional 16.0.0.82 (68) software.

Figure 1 Chain E of S-protein

Figure 2 Chain A of PL pro

2.2. Preparation of target proteins

The target proteins were prepared with the BIOVIA Discovery studio Visualizer v21.1.0.20298 software [20]. The Chain E (S-protein) and Chain A (PL pro) were used for the study and the other chains were deleted from the target proteins. The water molecules and heteroatoms were also removed from the target proteins. Now the polar hydrogen was added to both the target proteins. Then the active binding sites were defined with the ‘define and edit binding site’ feature of Discovery Studio Software and active site coordinates of S-protein (x=-35.41, y=8.55, z=29.13) and PL pro (x= 36.56, y= 9.35, z= 16.36) were saved for future use. The prepared proteins were saved in the PDB file format for future use.

2.2.        Preparation of compound library

The structure of the fourteen phytoconstituents of B. retusa was prepared manually with the Chem Draw Professional 16.0.0.82 (68) software. These constituents were already reported as potent against SARS-CoV-2 3CLPro protein [17]. Further, the SMILES ID of the structures was also saved for future use. The prepared compounds were saved in MDL SD File (*SDF) format. The structure of the standard drug Bepridil has been obtained from the PubChem database (PubChem CID 2351) and it was saved in the *SDF format.

2.3.        Molecular Docking Simulation study

The molecular docking simulation study (MDSS) was carried out with Autodock Vina on the virtual screening tool PyRx 0.8 software [21]. The prepared proteins were loaded in the 3D scene in the virtual platform of the software and converted to the PDBQT file format when made into macromolecules [22]. The originally downloaded proteins were reloaded and the unnecessary chains were removed from the scene except for Chain E (S-protein) and Chain A (PL pro). The sequence of the amino acids and the co-crystal ligand were revealed on expanding Chain E and Chain A of the respective proteins. The atoms of the co-crystal ligands were labeled to determine the accurate location of the co-crystal inhibitors present at the binding site of the protein. In the Vina search space of the PyRx 0.8 tool, the pre-defined active binding site coordinates were used to adjust the alignment of the 3D affinity grid box so that all the amino acids get covered at the active binding site of the protein. The size of the 3D affinity grid box was kept at default ay 25 Ã…. Finally, according to the standard protocols of the PyRx tool, MDSS were carried out [23].

2.4.        Visualisation and analysis of ligand interactions

The 2D interaction of the phytoconstituents with the best binding affinities with both the proteins was visualized with the Discovery Studio Visualizer Software. The 2D ligand interactions of the co-crystal-proteins was also visualized. The compound that did not form any conventional hydrogen bond with the active site residues was discarded from the study.

2.5.        Molecular Dynamics Simulations Studies

The molecular dynamics (MD) simulations studies was used to predict the most stable protein-ligand complex based on the vales of root mean square fluctuation (RMSF). A protein-ligand complex is considered to be stable if the RMSF value of all the amino acids is lower than 2.0 Ã… [24,25]. A protein-ligand complex of each constituent with S-protein and PL pro were generated and saved in the PDB file format. The MD simulations was carried out for only proteins along with the prepared complexes. The MD simulations was performed on the CABS Flex 2.0 server which utilizes the coarse-grained simulations of the protein motion [26]. All the parameters of the MD simulation study were kept at default.

3.            Results

3.1.        Features of the target proteins

The crystal structure of the chain E of ‘S’ protein and chain A of ‘PL Pro’ protein were retrieved from the RCSB-PDB website, shown in Figure 3 and Figure 4 respectively. The S-protein is made up of two chains, i.e, Chain A (sequence length 603), and Chain E (sequence length 229). The Chain E of the protein is complexed with co-crystal inhibitor viz. NAG (2-acetamido-2-deoxy-beta-D-glucopyranose). The PL pro is composed of one chain, i.e, Chain A (sequence length 315), complexed with HBA (p-hydroxy benzaldehyde).

Figure 3 Chain E of S-protein with co-crystal ligand

Figure 4 Chain A of papain-like protease with co-crystal ligand

3.2.       Molecular Docking Simulation study and visualization of ligand interactions

MDSS is an effective and competent tool for in-silico screening of drug compounds. A docking study on the PyRx tool can provide a binding affinity value (-kcal/mol) for each ligand so that the binding potential of ligands toward a protein can be ranked. The binding affinity of the compounds towards their respective binding site in the proteins is given in Table 1. For each ligand, the PyRx tool generates a maximum of 9 binding poses at the active binding site of the target protein.

Table 1 Binding affinity of each drug towards the active binding site of the S-protein and the PL pro

Drug	Binding affinity (-kcal/mol) S-protein	Binding affinity (-kcal/mol) PL pro
Co-crystal inhibitor	4.8	6.2
Bepridil (Standard)	6	7.2
(-)-Isochaminic acid	5.3	7.6
(+)- Sesamin	8.5	6.2
β-Sitosterol	7.7	6.7
Lupeol	7.7	6.2
(R)4-(1,5-dimethyl-3-oxo-4-hexenyl)-benzoic acid	6.4	8.2
(E)4-(1,5-dimethyl-3-oxo-1-hexenyl)-benzoic acid	6.5	7.8
(E)4-(1,5-dimethyl-3-oxo-1,4-hexadienyl)-benzoic acid	6.7	8.4
Elemicin	5.3	5.7
Stigmasterol	8.1	6.7
Gallic acid	4.7	6.7
Friedelin	8.6	6.6
Cumic acid	5.5	7.6
Ellagic acid	6.3	8.3
Ar-Todomatuic acid	6.1	7.8

Visualization of the 2D ligand interactions was done with Discovery Studio Visualizer software. The images of the 2D ligand interactions of the phytoconstituents with S-protein and PL pro are given in Figure 5 and Figure 6 respectively. The summary of the ligand interactions of each drug with the amino acids of both the proteins is also given in Table 2 and Table 3.

Figure 5 Visualization of 2D-interaction of S-protein with (A) NAG-(Co crystal), (B) Bepridil, (C) (+)-Sesamin, (D) (E)-4-(1,5-dimethyl-3-oxo-1,4-hexadienyl) benzoic acid, (E) (E)-4-(1,5-dimethyl-3-oxo-1-hexenyl) benzoic acid, (F) (R)-4-(1,5-dimethyl-3-oxo-4-hexenyl) benzoic acid, (G) Ar-todomatuic acid, (H) β-Sitosterol, (I) Ellagic acid, (J) Stigmasterol

Table 2 Summary of ligand interaction with the active site of S-protein

Figure 6 Visualization of 2D-interaction of PL Pro-protein with (A) HBA-(Co crystal), (B) Bepridil, (C) (-)-Isochaminic acid, (D) (R)-4-(1,5-dimethyl-3-oxo-4-hexenyl) benzoic acid, (E) (E)4-(1,5-dimethyl-3-oxo-1,4-hexadienyl)-benzoic acid, (F) (E)-4-(1,5-dimethyl-3-oxo-1-hexenyl) benzoic acid, (G) Cumic acid, (H) Ellagic acid, (I) Ar-todomatuic acid.

Table 3 Summary of ligand interaction with the active site of PL pro

Drugs	Conventional Hydrogen bond	Other interactive sites
HBA (Co-crystal inhibitor)	VAL11 (2.05Ã…)	VAL11(5.17Ã…), PRO59(5.40Ã…), TYR72 (4.86Ã…), LEU80 (4.09Ã…)
Bepridil	--	PRO77 (4.14Ã…), ALA68 (4.37Ã…, 4.40 Ã…), THR75 (3.65Ã…), ARG65 (4.47Ã…), PRO59 (4.50Ã…, 5.06Ã…, 3.44Ã…), LEU58 (4.36Ã…)
(-)-Isochaminic acid	--	PRO77 (4.46Ã…), PRO59 (5.29Ã…), ALA68 (4.48Ã…), PHE79 (5.09Ã…, 5.17Ã…), TYR72 (5.40Ã…, 4.07Ã…), LEU80 (4.21Ã…)
(R)4-(1,5-dimethyl-3-oxo-4-hexenyl)-benzoic acid	ARG65 (1.91Ã…)	ARG65 (4.84Ã…), PHE79 (5.15Ã…), PRO59 (3.98Ã…), LEU80 (4.10Ã…), ALA68 (4.25Ã…), TYR72 (4.21Ã…), THR75 (3.81Ã…)
(E)4-(1,5-dimethyl-3-oxo-1-hexenyl)-benzoic acid	--	PRO77 (3.73Ã…), PRO59(4.00Ã…), VAL11 (4.19Ã…), LEU80 (3.76Ã…)
(E)4-(1,5-dimethyl-3-oxo-1,4-hexadienyl)-benzoic acid	ARG65 (1.92Ã…), LEU80 (2.04Ã…)	ARG65 (4.41Ã…), LEU80 (4.19Ã…), PRO59 (4.41Ã…), ALA68 (4.47Ã…), THR75(3.74Ã…), VAL11 (4.97Ã…), TYR72 (4.26Ã…)
Cumic acid	ASP76 (1.98Ã…)	PHE79 (5.40Ã…), ALA68 (4.41Ã…), LEU80 (5.27Ã…, 5.03Ã…), PRO59 (4.20Ã…), TYR56 (5.20Ã…), TYR72 (4.19Ã…), TYR83 (4.72Ã…)
Ellagic acid	LEU80 (1.98Ã…), THR74 (2.45Ã…)	ALA68 (4.12Ã…, 4.20Ã…), THR75 (3.83Ã…), PRO59 (3.57Ã…, 4.81Ã…, 3.76Ã…, 3.93Ã…), PRO77 (5.23Ã…, 5.45Ã…)
Ar-Todomatuic acid	ARG65 (1.96Ã…, 2.02Ã…)	ARG65 (4.62Ã…), THR75 (3.72Ã…), ALA68 (4.31Ã…), PRO59 (4.42Ã…), PHE79 (4.89Ã…), LEU80 (4.33Ã…), TYR 72 (4.01Ã…)

3.3.       Molecular Dynamic Simulation Study

The MD simulation studies performed in the present study generated the RMSF values for each amino acid residue of the protein thereby furnishing an idea on the stability of each amino acid under a given set of parameters for a period of 10 nanoseconds. This study helps to validate the conformational stability of the protein-ligand complexes. The RMSF plot of the amino acids of S-protein and PL pro, both alone and in complex with all the phytoconstituents are given in Figure 7 and Figure 8 respectively.

4.            Discussion

MDSS is a computational technique for searching for an appropriate ligand that fits both energetically and geometrically the binding site of a protein [27]. A more negative binding affinity value suggests a better binding between a compound and a protein [23]. A low binding affinity value also indicates the low energy requirement for protein-ligand binding [27]. In all cases, the first pose is considered the best pose since it has the highest binding affinity and the last pose shows the lowest binding affinity towards the target protein.

Of all the compounds, 10 phytoconstituents of B. retusa have shown better binding affinity towards Chain E of the S-protein as compared to the binding affinities of the standards. These 10 constituents include (+)-Sesamin (-8.5 kcal/mol), (R)4-(1,5-dimethyl-3-oxo-4-hexenyl)-benzoic acid (-6.4 kcal/mol), (E)4-(1,5-dimethyl-3-oxo-1-hexenyl)-benzoic acid (-6.5 kcal/mol), (E)4-(1,5-dimethyl-3-oxo-1,4-hexadienyl)-benzoic acid (-6.7 kcal/mol), Stigmasterol (-8.1 kcal/mol), Friedelin (-8.6 kcal/mol), Ellagic acid (-6.3 kcal/mol), Ar-Todomatuic acid (-6.1 kcal/mol). Towards Chain A of PL pro, 7 phytoconstituents have shown better binding affinity than the standards. These include (-)-Isochaminic acid (-7.6 kcal/mol), (R)4-(1,5-dimethyl-3-oxo-4-hexenyl)-benzoic acid (-8.2 kcal/mol), (E)4-(1,5-dimethyl-3-oxo-1-hexenyl)-benzoic acid (-7.8 kcal/mol), (E)4-(1,5-dimethyl-3-oxo-1,4-hexadienyl)-benzoic acid (-8.4kcal/mol), Cumic acid (-7.6 kcal/mol), Ellagic acid (-8.3 kcal/mol), Ar-Todomatuic acid (-7.8 kcal/mol). These phytoconstituents are taken for visualization of ligand interactions and MD simulations study.

Figure 7 RMSF (Ã…) fluctuation plot of all the amino acid residues of Chain E of S-protein alone and in complex with all the phytoconstituents

Figure 8 RMSF (Ã…) fluctuation plot of all the amino acid residues of Chain A of PL pro alone and in complex with all the phytoconstituents

From the ligand interactions in Figures 5, 6, and Tables 2, and 3 above, it can be observed that for the S-protein, only four out of a total of sixteen compounds (including standards) has formed a conventional hydrogen bond with different amino acids at the active site. Bepridil, (+)- sesamin, (E)-4-(1,5-dimethyl-3-oxo-1,4-hexadienyl) benzoic acid, (E)-4-(1,5-dimethyl-3-oxo-1-hexenyl) benzoic acid, β-Sitosterol, stigmasterol were not able to form hydrogen bond. Therefore, the remaining four compounds were taken for further analysis. Similarly, for the PL pro, most of the phytoconstituents have formed a conventional hydrogen bond with different amino acids of the active site. The three compounds namely, Bepridil, (-)-isochaminic acid, and (E)4-(1,5-dimethyl-3-oxo-1-hexenyl)-benzoic acid were not able to form any hydrogen bond. Therefore, the remaining six phytoconstituents were taken for the MD simulations study.

For each amino acid, a lower RMSF value indicates limited flexibility and a higher RMSF value indicates high flexibility in a given system [25]. For the S-protein, the RMSF value for 78.75% of all the amino acids of the protein without the presence of any ligand was found to be below 2 Ã…. For the S-NAG complex, the RMSF value of a total of 78.75% amino acids was found to be lower than 2 Ã…. For the S- (R)-4-(1,5-dimethyl-3-oxo-4-hexenyl) benzoic acid complex, S- Ar-todomatuic acid complex, S- Ellagic acid complex, 78.75% amino acids have exhibited RMSF value lower than 2 Ã…. For the PL Pro protein, without the presence of any ligand, the RMSF value of a total of 92.18% amino acids was found to be lower than 2 Ã…. For PL Pro- HBD complex, the RMSF value of a total of 92.51% amino acids was found to be lower than 2 Ã…. For PL Pro- Ellagic acid complex, the RMSF value of a total of 93.16% amino acids was found to be lower than 2 Ã…. For PL Pro- Cumic acid, the RMSF value of a total of 94.46% amino acids was found to be lower than 2 Ã…. For PL Pro- Ar-todomatuic acid complex, the RMSF value of a total of 93.49% amino acids was found to be lower than 2 Ã….

For PL Pro- (R)4-(1,5-dimethyl-3-oxo-4-hexenyl)-benzoic acid complex, the RMSF value of a total of 97.07% amino acids was found to be lower than 2 Ã…. For PL Pro- (E)4-(1,5-dimethyl-3-oxo-1,4-hexadienyl)-benzoic acid complex, the RMSF value of a total of 94.46% amino acids was found to be lower than 2 Ã…. The constituents of B. retusa namely, (R)4-(1,5-dimethyl-3-oxo-4-hexenyl)-benzoic acid, Ellagic acid, Ar-Todomatuic acid were able to form the same stable complexes similar to that of its co-crystal inhibitors of the S-protein. Phytoconstituents such as cumic acid, E-4-(1,5-dimethyl-3-oxo-1,4-hexadienyl)-benzoic acid has been able to form a more stable complex with the PL pro than that of its co-crystal inhibitor. Of these (R)4-(1,5-dimethyl-3-oxo-4-hexenyl)-benzoic acid is found to form the most stable complex with PL pro.

Nowadays, researchers have used new research methods that deviates from the traditional herbal research. For example, in-silico techniques such as MDSS and MD simulations have been increasingly applied in drug discovery research to identify promising phytocompounds for the treatment of various diseases [17, 28-39]. Medicinal plants are an important source of clinically important phytocompounds and many pharmaceutical drugs have been developed from traditional herbal remedies [40-44]. However, some phytocompounds have poor oral bioavailability. To overcome the bioavailability issues associated with natural products, novel drug delivery systems have been adopted by many researchers as the solution [45, 46]. Artificial intelligence, machine learning, and supervised machine learning have also been utilized in the process of drug discovery and development [47, 48]. With the advancement in science, pharmaceutical researchers are using new techniques in drug discovery programs. In the present study, we have also utilized in-silico techniques that is sustainable, safe, and economical for the identification of promising bioactive molecules against SARS-CoV-2.

5.            Conclusions

The present in-silico study has revealed that the phytoconstituent (R)4-(1,5-dimethyl-3-oxo-4-hexenyl)-benzoic acid of B. retusa is a potential inhibitor of the S-protein and PL pro of the SARS-CoV-2. It has exhibited better binding affinity towards the active binding site of the S-protein and PL pro in comparison to the standards. (R)4-(1,5-dimethyl-3-oxo-4-hexenyl)-benzoic acid has shown similar molecular interactions with that of the standards in the active binding sites of both the proteins. However, the present study is confined to only the in-silico model, hence, further, (in-vitro/in-vivo) study is needed to determine the complete inhibitory potential of the phytocompound against S-protein and PL pro of SARS-CoV-2.