# Structure-Based Virtual Screening, Docking, ADMET, Molecular Dynamics, and MM-PBSA Calculations for the Discovery of Potential Natural SARS-CoV-2 Helicase Inhibitors from the Traditional Chinese Medicine

Authors: Ahmed M. Metwaly; Alaa Elwan; Abdul-Aziz M. M. El-Attar; Sara T. Al-Rashood; Ibrahim H. Eissa
Journal: Journal of Chemistry
 (2022)

---

## Abstract

Continuing our antecedent work against COVID-19, a set of 5956 compounds of traditional Chinese medicine have been virtually screened for their potential against SARS-CoV-2 helicase (PDB ID: 5RMM). Initially, a fingerprint study withVXG, the ligand of the target enzyme, disclosed the similarity of 187 compounds. Then, a molecular similarity study declared the most similar 40 compounds. Subsequently, molecular docking studies were carried out to examine the binding modes and energies. Then, the most appropriate 26 compounds were subjected to in silico ADMET and toxicity studies to select the most convenient inhibitors to be: (1R,2S)-ephedrine (57), (1R,2S)-norephedrine (59), 2-(4-(pyrrolidin-1-yl)phenyl)acetic acid (84), 1-phenylpropane-1,2-dione (195), 2-methoxycinnamic acid (246), 2-methoxybenzoic acid (364), (R)-2-((R)-5-oxopyrrolidin-3-yl)-2-phenylacetic acid (405), (Z)-6-(3-hydroxy-4-methoxystyryl)-4-methoxy-2H-pyran-2-one (533), 8-chloro-2-(2-phenylethyl)-5,6,7-trihydroxy-5,6,7,8-tetrahydrochromone (637), 3-((1R,2S)-2-(dimethylamino)-1-hydroxypropyl)phenol (818), (R)-2-ethyl-4-(1-hydroxy-2-(methylamino)ethyl)phenol (5159), and (R)-2-((1S,2S,5S)-2-benzyl-5-hydroxy-4-methylcyclohex-3-en-1-yl)propane-1,2-diol (5168). Among the selected 12 compounds, the metabolites, compound 533 showed the best docking scores. Interestingly, the MD simulation studies for compound 533, the one with the highest docking score, over 100 ns showed its correct binding to SARS-CoV-2 helicase with low energy and optimum dynamics. Finally, MM-PBSA studies showed that 533 bonded favorably to SARS-CoV-2 helicase with a free energy value of −83 kJ/mol. Further, the free energy decomposition study determined the essential amino acid residues that contributed favorably to the binding process. The obtained results give a huge hope to find a cure for COVID-19 through further in vitro and in vivo studies for the selected compounds.


---

## Body

## 1. Introduction

The WHO disclosed on December 25, 2021 that the confirmed COVID-19 infections globally became over 276 million including more than 5 million dead persons [1]. These massive numbers demand enormous work from scientists all over the world to find a cure.The utilization of natural products for the treatment has been mentioned since the oldest historical points [2]; the traditional medicines were unlimited sources for bioactive natural compounds such as flavonoids [3–5], alkaloids [6], saponins [7–9], isochromenes [10], α-pyrones [11, 12], diterpenoids [13], sesquiterpenoids [14, 15], and steroids [16].Traditional Chinese medicine (TCM) is an ethnomedicine that authenticates the experience of ancient Chinese people in the treatment of different illnesses [17]. TCM is a reflection of a great experience of clinical practice that extended for thousands of years [18]. To date, the TCM remedies are still utilized effectively in China as well as several places of the world [19, 20].Computer-aided drug design methodologies play an ever-increasing essential role in the discovery of new drugs [21, 22]. These methodologies have been very effective in the identification of new promising drug candidates with a noticeable limitation in time, cost, effort, and use of animal models [13–15]. The application of in silico methodologies included molecular docking [16–23], molecular design [24, 25], rational drug design [26–31], computational chemistry [32, 33], toxicity [34–36], ADMET [37–39], and DFT [40] assessments.Our teamwork utilized computer-based methodologies to determine potential inhibitors against COVID-19 in various reports. For example, metabolites ofArtemisia sublessingiana [41]‏ and Monanchora species [42] were examined in silico against COVID-19. We suggested four isoflavonoids between a set of 59 as the most promising inhibitors against hACE2 and Mpro [43]. Recently, our team adjusted a multistep in silico filtration technique to select the most promising compound through a huge group of compounds against a certain COVID-19 protein. For instance, vidarabine was selected to be the most promising inhibitor against SARS-CoV-2 nsp10 [44]. Complementarily, the most convenient semisynthetic molecule against PLpro has been determined via a set of 69 compounds [45].Helicases are pivotal enzymes in the viral lifecycle because of their responsibility to separate the dsDNA or RNA strands as well as their essential role in the process of RNA replication and repair [46]. Helicases can translocate molecules along the double-stranded (ds) DNA as well as RNA in a certain direction. Additionally, it can unwind (separate) the complementary strand of the DNA duplex through the dissociation of the hydrogen bonds between the nucleotide bases [47].In this work, a collection of 5956 natural compounds, that were derived from traditional Chinese medicine and available athttp://tcm.cmu.edu.tw/, has been subjected to structure and ligand-based in silico approaches (Figure 1) to determine the most convenient SARS-CoV-2 helicase inhibitors. The starting step in our research was (3S,4R)-1-acetyl-4-phenylpyrrolidine-3-carboxylic acid (VXG), the co-crystallized ligand of the SARS-CoV-2 helicase (PDB ID: 5RMM). VXG showed a high binding affinity against the target enzyme. Accordingly, it is expected according to the SAR principles that any compound with a similar structure could have a high binding affinity too. The utilized in silico methods included molecular structure similarity and fingerprint study against the VXG. Then, molecular docking against SARS-CoV-2 helicase (PDB ID: 5RMM) was conducted to examine the binding. ADMET and toxicity were utilized to make sure about the likeness of the selected compounds. Finally, molecular dynamics (MD) simulation experiments (RMSD, RMSF, Rg, SASA, and H-bonding) over 100 ns for the compound of the highest docking score, as well as MM-PBSA studies, were preceded to confirm the correct binding mode.Figure 1
The employed computational techniques.


## 2. Results and Discussion

### 2.1. Molecular Filtration Using Fingerprint Method

The ligand-basedin silico approach depends on the computation of chemical and physical properties of a molecule (ligand) and comparison of these properties with some biologically active compounds [48]. The fingerprint study is one of the ligand-based in silico methods that is vastly employed to predict the chemical structure’s similarity or dissimilarity of two compounds or more [49, 50]. During the fingerprint study, the computer converts the chemical descriptors of a molecule to mathematical symbols. The obtained data are presented as bit strings. These strings describe the presence (1) or absence (0) of a certain 2D fragment or atomic descriptor (property) in the test and reference molecules [51, 52]. The co-crystallized ligand is a molecule that has a very high affinity to bind to a specific protein forming a ligand-protein complex in a crystallized form [53]. In consequence, the chemical structure of the co-crystallized ligand could be utilized effectively as a starting point to design and discover a potential inhibitor against the target protein.Discovery studio 4.0 software was employed to examine the fingerprint similarity of 5956 natural compounds, which were derived from traditional Chinese medicine, againstVXG, the co-crystallized ligand of SARS-CoV-2 helicase (PDB ID: 5RMM). The experiment determined 187 compounds to be the most similar candidates to VXG (Table 1). The study compared the following descriptors (properties) in the chemical structures of the experiment set and VXG: H-bond acceptor [54], H-bond donor [55], charge [56], hybridization [57], positive ionizable atoms [58], negative ionizable atoms [59], halogens [60], aromatic groups [61], and aligning with the ALogP [62] of fragments as well as atoms.Table 1
Fingerprint similarity between 187 natural compounds andVXG.
Comp.SimilaritySASBSCComp.SimilaritySASBSCVXG11660054120.582205602580−1481381160100615850.581927586553−1321671169118−315660.581457439301152151180166−148100.581121197173−3125811531281354490.5808634312882334711431012323790.580247329113125364190−77651540.579775258−919637911431012332390.579154589563−13538011431012323560.57854523450−694111112425432580.57820545132634451186146−2035890.578089496404−424961196201−303620.577495272171825011163145354390.577398608599−1545071158136832550.577267452329252611431012315840.576874454333053311337033710.576792169127−35521173123−76230.576792169127−35541173121−7380.57608726561895551173121−754410.576077602591−1485771154901219820.575087494405−40610111639503980.57485961706191156131101140.5611111011465637116210848160.5611234547331169157−3920.556886931737561206193−401690.556886931737821220219−5451570.55670110828587921207178−414510.554878182162−167941196177−3051540.554307148101188031169147−37680.55418179157−138061210206−441900.55329910931578071183144−175570.55298167136−18181999676290.55176154−1023791213229−474950.54333316313434050.7771369304330.538462913754420.7381276397540.537764178165−12850.7311255415660.534743177165−112600.722117−4493540.534535178167−122080.69712412425970.53271171155−5910.68912211445960.53114816213941020.68912211446010.528302168152−21950.681109−6577520.527473192198−262800.6741166505680.526946176168−103440.65812118457900.526012182180−163700.658121184523680.523041227268−616720.65216790−14390.522222188194−22480.6491138531580.521978951671580.64911385351520.521898143108231000.646104−5622500.520408153128132460.644105−3615910.5197371581388840.638118194851530.517699117604951690.6361598471500.517241908768170.63312633402970.515789147119193960.62410646023700.515738213247−472450.622102−26451670.51567181185−15570.621065604130.51398614712019650.621065604910.5136991501261619420.614339497355−4351770.512195189203−2354150.614108592510−1384150.5119451501271654340.614017587502−13351510.511696175176−939410.6138334262402823670.511364225274−5940500.613551489343−355440.511299181188−1554410.613124626567−17251590.51111192147450390.613119458293−42660.510903164155240550.612984491347−373740.5105631451182135910.612751519393−654180.5104171471221923480.612622563465−10951550.510288124774239120.612319507374−532070.5102741491261735790.612128535420−8151490.51005203232−3739190.6114294282462651800.508197186200−2050550.611111440266142990.5070421441182251680.6115893851580.50632912071465390.609162100451760.505208194218−283420.60811015563970.505051150131165110.608178127−126140.505017151133157840.607193152−276220.504983152135145640.605182135−165120.50461516415922010.604186142−2051730.50358211253−455420.589905187151−216390.503356150132164920.58885169121−3500.50316515915072280.586319180141−141160.502982253337−8723650.584081521438−67800.502092120734632540.58408587551−13332880.502092240312−742060.583178624616−1704540.502075121754515710.582915464342−104190.5017421441212235910.582851503409−496210.501651521371415930.582512473358−192030.501449173179−7590.5823539946732910.501006249331−83640.58235399467860.515414212860.582278276201781180.51003466SA, bits number in both similar natural compounds and VXG; SB, bits number in the similar natural compounds but not in VXG; SC, bits number in VXG but not in the similar natural compounds.


### 2.2. Molecular Similarity

The difference between molecular similarity and fingerprint studies is that the fingerprint study computes the presence and/or absence of specific 2D atom paths and descriptors regarding fragments or substructures in the examined compounds [63]. Contrastingly, the molecular similarity calculates specific molecular descriptors considering the whole chemical structure of compounds. These descriptors are steric, topological, electronic, and/or physical [64]. Utilizing Discovery studio 4.0 software, a molecular similarity study was done on the most similar 187 compounds against VXG. The applied descriptors in this study (Figure 2 and Table 2) were partition coefficient (ALogp) [65], molecular weight (M. W) [66], H-bond donors (HBA) [67], H-bond acceptors (HBD) [68], rotatable bond numbers [69], number of rings as well as aromatic rings [70], and minimum distance [71] together with the molecular fractional polar surface area (MFPSA) [72]. The study was adapted to select the most similar 40 metabolites (Figure 3).Figure 2
Molecular similarity of the examined compounds andVXG.Table 2
Molecular descriptors of the examined 40 compounds andVXG.
Comp.ALogpM. WtHBAHBDRotatable bondsRingsAromatic ringsMFPSAMinimum distanceVXG0.71233.26312210.2374050.62219.24323210.3070.357842.12205.25313210.1870.4468162.1221.3322210.1810.4544331.69178.19313110.2390.5291001.09149.19212110.2520.5313422.76204.22312210.2240.5415331.97274.27514210.2290.5528181.53195.26323110.1870.5533641.44152.15312110.2850.5552461.91178.19313110.240.559590.8151.21222110.2640.56640.8151.21222110.2640.565391.57286.32423310.2390.5641501.54178.19322110.30.568852.03191.23201210.1430.5753981.81164.2213110.2020.5961951.45148.16202110.210.597571.23165.23223110.1650.614651.23165.23223110.1650.6148172.34205.3212210.1020.64851682.17276.37334210.1980.67251532.59249.35333210.1810.67851591.56195.26334110.2310.6911182.22203.28200210.1220.6953742.4317.36422310.2440.7062451.93162.19203110.1430.713481.77179.26213110.1050.715581.77179.26213110.1050.7151142.2194.23304110.1570.7316103.51230.26312210.1920.7383882.01333.36523310.2640.74351693.05260.37223210.1360.762601.93164.2204110.1410.765912.99188.22201210.1310.7691022.99188.22201210.1310.7692801.98177.24201210.0620.77551553.04263.38333210.1710.7826371.37336.77533310.2820.7842082.16191.27201210.0570.795800.9292.28631310.3370.81Figure 3
The most similar compounds toVXG.


### 2.3. Docking Studies

Computer-aided drug design applies various techniques to optimize natural products into potentially active against certain biological targets [44, 45]. A molecular docking technique was applied for 40 VXG similar compounds against SARS-CoV-2 helicase (PDB ID: 5RMM). The binding modes and affinities of these compounds were examined.The target protein was downloaded from Protein Data Bank (http://www.pdb.org), and Molecular Operating Environment (MOE .14) was used for the docking analysis. The docking process was validated by the re-docking of VXG inside the active pocket of the helicase protein. The root mean square deviation (RMSD) between the re-docked and the co-crystallized conformers was 0.74 A, which confirms the validity of the docking protocol (Figure 4).Figure 4
(a) 3D and (b) 2D superimposition of the re-docked conformer ofVXG over the co-crystallized one with an RMSD value of 0.74 A.
(a)(b)Out of the examined 40 compounds, 26 displayed correct binding modes as well as good free energies. The promising compounds are as follows: (+)-methylpseudoephedrine (48), (1R,2S)-ephedrine (57), (1R,2S)-N-methylephedrine (58), (1R,2S)-norephedrine (59), (3S)-2,2-dimethyl-3,5-dihydroxy-8-hydroxymethyl-3,4-dihydro-2H,6H-benzo[1,2-b5,4-b′]dipyran-6-one (80), 2-(4-(pyrrolidin-1-yl)phenyl)acetic acid (84), (4S,5R)-ephedroxane (85), 1-phenylpropane-1,2-dione (195), 2-methoxycinnamaldehyde (245), 2-methoxycinnamic acid (246), 2-methoxybenzoic acid (364), 3′-o-acetylhamaudol (374), 3′-o-propionylhamaudol (388), (R)-2-((R)-5-oxopyrrolidin-3-yl)-2-phenylacetic acid (405), (Z)-6-(3-hydroxy-4-methoxystyryl)-4-methoxy-2H-pyran-2-one (533), 6,7-dihydroxy-2-(2-phenylethyl)-5,6,7,8-tetrahydrochromone (539), 7-demethylsuberosin (610), 8-chloro-2-(2-phenylethyl)-5,6,7-trihydroxy-5,6,7,8-tetrahydrochromone (637), 3-((R)-hydroxy ((S)-1-methylpiperidin-2-yl)methyl)phenol (816), (R)-((S)-1-methylpiperidin-2-yl) (phenyl)methanol (817), 3-((1R,2S)-2-(dimethylamino)-1-hydroxypropyl)phenol (818), (R)-4-(1-Hydroxy-2-(methylamino)ethyl)-7,7-dimethyl-5,6,7,8-tetrahydronaphthalen-1-ol (5153),(R)-4-(1-hydroxy-2-(methylamino)ethyl)-8,8-dimethyl-6,7,8,9-tetrahydro-5H-benzo [7]annulen-1-ol (5155), (R)-2-ethyl-4-(1-hydroxy-2-(methylamino)ethyl)phenol (5159), (R)-2-((1S,2S,5S)-2-benzyl-5-hydroxy-4-methylcyclohex-3-en-1-yl)propane-1,2-diol (5168), and (1S,4R,5S)-4-benzyl-5-(2-hydroxypropan-2-yl)-2-methylcyclohex-2-en-1-ol (5169).The docking scores of the experienced ligands are predicted and summarized in Table3. The binding modes of the tested ligands inside the active site of the target protein were depicted. The binding poses of the top five compounds with the highest energy scores as well as the most perfect modes were selected for detailed discussion as representative examples.Table 3
∆G (in kcal/mole) of the most similar 40 compounds toVXG.
CompoundName∆G (kcal/mole)48(+)-Methylpseudoephedrine−15.9257(1R,2S)-Ephedrine−15.6458(1R,2S)-N-methylephedrine−15.2459(1R,2S)-Norephedrine−14.2964(1S,2S)-Norpseudoephedrine−13.6565(1S,2S)-Pseudoephedrine−11.6680(3S)-2,2-Dimethyl-3,5-dihydroxy-8-hydroxymethyl-3,4-dihydro-2H,6H-benzo[1,2-b5,4-b′]dipyran-6-one−14.27842-(4-(Pyrrolidin-1-yl) phenyl) acetic acid−15.4685(4S,5R)-Ephedroxane−14.0191(E)-3-Butylidene phthalide−13.31100(S)-Cathinone−12.6102(Z)-3-Butylidene phthalide−13.91141-(2,4-Dimethoxyphenyl)-1-propanone−13.531181,3,6,6-Tetramethyl-6,7-dihydroisoquinolin-8(5H)-one−12.41150(4S,5R)-Ephedroxane−13.021951-Phenylpropane-1,2-dione−11.072082,3,4-Trimethyl-5-phenyloxazolidine−13.762452-Methoxycinnamaldehyde−14.922462-Methoxycinnamic acid−14.222602-Phenylethyl acetate−13.292803,4-Dimethyl-5-phenyloxazolidine−13.213423-Butylidene-4-hydro-phthalide−13.243642-Methoxybenzoic acid−12.163743′-o-Acetylhamaudol−18.213883′-o-Propionylhamaudol−18.133984-(4-Hydroxyphenyl)-2-butanone−12.88405(R)-2-((R)-5-Oxopyrrolidin-3-yl)-2-phenylacetic acid−15.744334-Hydroxy-3-methoxycinnamaldehyde−13.11533(Z)-6-(3-Hydroxy-4-methoxystyryl)-4-methoxy-2H-pyran-2-one−17.15396,7-Dihydroxy-2-(2-phenylethyl)-5,6,7,8-tetrahydrochromone−18.366107-Demethylsuberosin−15.486378-Chloro-2-(2-phenylethyl)-5,6,7-trihydroxy-5,6,7,8-tetrahydrochromone−18.858163-((R)-Hydroxy ((S)-1-methylpiperidin-2-yl)methyl)phenol−17.09817(R)-((S)-1-Methylpiperidin-2-yl) (phenyl)methanol−15.958183-((1R,2S)-2-(Dimethylamino)-1-hydroxypropyl)phenol−15.795153(R)-4-(1-Hydroxy-2-(methylamino)ethyl)-7,7-dimethyl-5,6,7,8-tetrahydronaphthalen-1-ol−17.845155(R)-4-(1-Hydroxy-2-(methylamino)ethyl)-8,8-dimethyl-6,7,8,9-tetrahydro-5H-benzo[7]annulen-1-ol−16.595159(R)-2-Ethyl-4-(1-hydroxy-2-(methylamino)ethyl)phenol−17.895168(R)-2-((1S,2S,5S)-2-Benzyl-5-hydroxy-4-methylcyclohex-3-en-1-yl)propane-1,2-diol−18.45169(1S,4R,5S)-4-Benzyl-5-(2-hydroxypropan-2-yl)-2-methylcyclohex-2-en-1-ol−17.42VXG(3S,4R)-1-Acetyl-4-phenylpyrrolidine-3-carboxylic acid−19.37Starting with the binding interactions and orientation of the co-crystallized ligand (VXG) inside the active SARS-CoV-2 helicase (PDB ID: 5RMM), it revealed a binding affinity value of −19.37 kcal/mol. It showed a characteristic four hydrogen bonding interactions through carboxylate moiety of pyrrolidine ring with the essential amino acids SER486, ASN516, and ASN177. In addition, two hydrophobic interactions were formed between pyrrolidine ring and amino acid residues HIS554 and TYR515 (Figure 5).Figure 5
(a) 3D and (b) 2D ofVXG docked into the active site of SARS-CoV-2 helicase.
(a)(b)The results of docking studies showed that the tested ligands have orientations and binding interactions similar to that ofVXG against SARS-CoV-2 helicase. Three- and two-dimensional representations of binding modes of the most potent derivatives 533, 637, 84, 195, and 364 inside the active site of the target protein are depicted in Figures 6–10.Figure 6
(a) 3D and (b) 2D images of the docked compound533 into the active site of SARS-CoV-2 helicase.
(a)(b)Figure 7
(a) 3D and (b) 2D images of the docked compound637 into the active site of SARS-CoV-2 helicase.
(a)(b)Figure 8
(a) 3D and (b) 2D images of the docked compound84 into the active site of SARS-CoV-2 helicase.
(a)(b)Figure 9
(a) 3D and (b) 2D images of the docked compound195 into the active site of SARS-CoV-2 helicase.
(a)(b)Figure 10
(a) 3D and (b) 2D images of docked compound364 into the active site of SARS-CoV-2 helicase.
(a)(b)Compound533 exhibited an interesting binding mode similar to that of the co-crystallized ligand against SARS-CoV-2 helicase with a docking score of −17.10 kcal/mol. It keeps the hydrogen bonding interactions with the essential amino acids SER486, ASN516, and ASN177. Also, it formed two additional hydrogen bonds with ASN179 and SER485 residues. Furthermore, compound 533 was incorporated in two hydrophobic interactions with amino acid residues HIS554 and TYR515 (Figure 6).For compound637, the binding affinity was −18.85 kcal/mol. Such compound exhibited the best binding mode into the target protein, where it completely occupied the protein through seven hydrogen bonding interactions with SER486, ASN516, ASN177, SER485, and HIS554 residues. In addition, the terminal phenyl ring formed three hydrophobic interactions with HIS554, LEU412, and ARG560 residues (Figure 7).Concerning the binding mode of compound84 against SARS-CoV-2 helicase, the binding energy was −15.46 kcal/mol. The structural similarity of that compound with the co-crystallized ligand revealed the same binding mode against the receptor, where three hydrogen bonds were formed with SER486, ASN516, and ASN177 and three hydrophobic interactions were molded with HIS554 and PRO408 (Figure 8).The binding affinity of compound195 was −11.07 kcal/mol. Such affinity was represented by four hydrogen bonds with the key amino acids SER486, ASN516, and ASN177 and one hydrophobic interaction with HIS554 (Figure 9).Finally, analyzing the binding interactions of compound364 indicated a binding score of −12.16 kcal/mol. The carboxylate moiety formed three hydrogen bonding interactions with the key amino acids SER486, ASN516, and ASN177 while the phenyl ring was incorporated in hydrophobic interaction with HIS55 (Figure 10).


### 2.4. ADMET Studies

The likeness of any molecule to be approved as a drug depends greatly on its pharmacokinetic properties as well as its activity. Subsequently, the investigation of the ADMET profile of a molecule should be considered in the early stages of drug design and discovery to avoid the withdrawal possibility of the drug from the pharmaceutical market [73]. These descriptors identify the absorption, distribution, metabolism, excretion, as well as the toxicity of the examined compound. Although, there are different in vitro experiments that can determine the ADMET profile, in silico determination is an available and reliable tool with the profit of being faster, cheaper, as well as and lifesaver of the experimental animals [74].The predicted ADMET profiles of the 26 compounds that showed correct modes of binding besides Remdesivir, the reference drug, are shown in Table4 and Figure 11. Compounds 48, 58, 85, 245, 610, 816, 817, 5153, 5155, and 5169 were excluded because of their predicted strong ability to pass the blood-brain barriers which may be combined with a CNS toxicity. Among the excluded compounds, compounds 48, 58, 816, 817, 5153, and 5155 were predicted to be inhibitors against the CYP2D6 enzyme which would cause hepatotoxicity. The predicted intestinal absorption and aqueous solubility of all compounds were good to optimal.Table 4
ADMET profile of the 26 compounds with the best docking scores.
Comp.BBB levelaHIAbAqcCYP2D6dPPBe48104TF57204FF58104TF59304FF80303FF84203FT85103FF195203FT245103FT246203FT364204FT374203FF388303FF405304FF533303FT539303FT610103FT637303FT816203TF817103TF818204FF5153203TT5155203TT5159304FF5168204FT5169103FTSimeprevir432FTaBBB, ability to pass the blood-brain barrier, 1 is high, 2 is medium, 3 is low, and 4 is very low; bHIA, human intestinal absorption level, 0 is good, 1 is moderate, 2 is poor, and 3 is very poor; cAq, aqueous solubility level, 0 is extremely low, 1 is very low, 2 is low, 3 is good, and 4 is optimal; dCYP2D6, inhibition of CYP2D6 enzyme, T is an inhibitor and F is a noninhibitor; ePPB, F means less than 90% and T means more than 90%.Figure 11
Results of the ADMET study.


### 2.5. Toxicity Studies

The early prediction of toxicity is a crucial step that minimizes drug failure because of toxicity in the development stage or the clinical trials [75]. In silico prediction of toxicity is a credible approach that plays an essential role in drug design and discovery of lead compounds because in vitro and in vivo approaches are usually controlled by strict ethical regulations, time, and availability of resources [76], whereas the in silico prediction is based on a structure-activity relationship toxicology. The software compares the essential structural descriptors of the examined compounds with a huge library of hundreds of thousands of reported safe and toxic compounds [77] (Supporting data (available here)). Discovery studio 4.0 software was employed to predict the toxicity profile of the selected compounds after the ADMET study against 7 models. The applied models are FDA rat carcinogenicity [78, 79], mouse carcinogenic potency (TD50) [80], rat maximum tolerated dose (MTD) [81, 82], rat oral LD50 [83], rat chronic LOAEL [84, 85], ocular irritancy, and skin irritancy [86]. According to the obtained results (Table 5), compounds 80, 388, and 539 were eliminated due to the predicted high carcinogenic potency.Table 5
In silico toxicity profile of 15 compounds with good ADMET profile.
Comp.FDA rat carcinogenicity (mouse-female)TD50 (mg·kg−1·day−1)MTD (g·kg−1)Rat oral LD50 (g·kg−1)LOAEL (g·kg−1)Ocular irritancySkin irritancySimeprevirNot a carcinogen2.01380.0029670.2088350.0021057MildNone57Not a carcinogen306.6850.1269250.6785570.1091SevereNone59Not a carcinogen209.0130.1306691.178220.153373SevereNone80Multicarcinogen30.96540.3205630.1126450.0104398ModerateNone84Not a carcinogen86.27450.4320430.6515370.0466934SevereNone195Not a carcinogen734.3760.09202340.803940.58711MildNone246Not a carcinogen896.4370.1859080.9757830.0690491MildMild364Not a carcinogen1,152.330.1781251.100170.269177MildNone388Multicarcinogen77.44010.1582120.1285760.0111578SevereMild405Not a carcinogen1,019.870.3558430.5655660.0958513ModerateNone533Not a carcinogen587.5160.1090020.7653060.026065MildMild539Multicarcinogen71.15820.123540.3662150.0229232SevereNone637Not a carcinogen30.33650.1365950.4284340.0201137SevereMild818Not a carcinogen145.390.3557360.712450.0937066SevereNone5159Not a carcinogen227.5990.5599770.5267870.0808567SevereNone5168Not a carcinogen128.9110.2428542.096460.0494677ModerateMild


### 2.6. Molecular Dynamics (MD) Simulations

Despite the ability of molecular docking studies to expect the mode of binding of a compound inside a specific protein correctly, it has a serious defect that it deals with the protein as a rigid unit. Resultantly, it does not compute the conformational changes that happen in the protein because of ligand binding [87]. On the contrary, the MD simulations can adequately describe the behavior of a protein at the atomic level in full detail and at very accurate temporal resolution [88]. In accordance, the MD simulations have the advantage of being able to predict the conformational changes that occurred in the protein after ligand binding [89]. Furthermore, MD simulation studies can effectively compute various factors related to the energy of the protein-ligand complex for a determined time. Subsequently, it accurately describes the binding mode, stability, and flexibility of the ligand inside the target protein [90].The first successful MD simulation experiment of a protein (bovine pancreatic trypsin inhibitor) was published in Nature in 1977 [91]. Fortunately, because of the recently introduced supercomputer hardware, especially the advanced graphics processing units, MD simulation experiments became much more accessible, powerful, and accurate [92].In the MD simulation study, the forces on every atom of the examined ligand-protein complex are computed at every ultra-short time interval according to the basics of the “force field” [93]. The computed force field can be utilized to describe the position and velocity of atoms at each time interval. The force field is a physical expression that describes the functional potential energy of atoms. The force field is calculated based on Newton’s laws of motion considering bonded interactions (bonds, angles, and dihedrals) in addition to nonbonded interactions (van der Waals potentials and Coulomb potentials) between all atoms of the complex. This step is repeated billions of times to produce the atomic trajectories for a specific time interval [94].Several MD simulation studies were employed to investigate the stability and mimic the dynamic of compound533, (Z)-6-(3-hydroxy-4-methoxystyryl)-4-methoxy-2H-pyran-2-one, that exhibited the best docking score inside SARS-CoV-2 helicase for 100 ns.First, the interaction of a ligand inside the active site of a protein leads to some changes in the structure of that protein [95]. In consequence, the dynamics and the conformational changes of the SARS-CoV-2 helicase-533 complex were computed as root mean square deviation (RMSD) to detect the stability due to binding. It is observed that the SARS-CoV-2 helicase and 533 exhibited lower RMSD with no major fluctuations indicating their greater stability (Figure 12(a)). Interestingly, the SARS-CoV-2 helicase-533 complex was stable till 90 ns∼. Although the SARS-CoV-2 helicase-533 complex showed a minor fluctuation later, it reached equilibrium again.Figure 12
M D simulation results: (a) RMSD values of SARS-CoV-2 helicase,533, and SARS-CoV-2 helicase-533 complex; (b) RMSF for SARS-CoV-2 helicase; (c) Rg of SARS-CoV-2 helicase; (d) SASA of SARS-CoV-2 helicase; and (e) H-bonding between SARS-CoV-2 helicase-533 complex in the MD run.
(a)(b)(c)(d)(e)Second, the flexibility of SARS-CoV-2 helicase was calculated in terms of root mean square fluctuation (RMSF) to calculate the differences in flexibility in the SARS-CoV-2 helicase-533 complex during the 100 ns of the MD simulations. The decrease of RMSF values during the MD simulation in 50–100 residue areas (Figure 12(b)) denotes that SARS-CoV-2 residues were more rigid and stabilized after binding to 533.Third, the radius of gyration (Rg) is a crucial parameter that is linked to the protein stability according to the change in its volume. Rg is defined as the root mean square distance (RMSD) of a weighted mass group of atoms from their mass center [96, 97]. Thus, the calculation of Rg identifies the dimensions as well as the compactness of the SARS-CoV-2 helicase-533 complex. The lower degree of fluctuation throughout the simulation period indicates the greater compactness of a system. The Rg of the SARS-CoV-2 helicase-533 complex was found to be lower than the starting period (Figure 12(c)) displaying compactness and stability.Fourth, the interaction between protein-ligand complexes and solvents was measured by solvent accessible surface area (SASA) over the simulation period. So, the SASA of the SARS-CoV-2 helicase-533 complex was calculated to provide the extent of the conformational changes that occurred during binding. Interestingly, SARS-CoV-2 helicase featured a reduction of the surface area showing a relatively lower SASA value than the starting period (Figure 12(d)).Finally, hydrogen bonding between a SARS-CoV-2 helicase-533 complex is essential to stabilize the structure. MD simulation studies showed that the highest number of conformations of the SARS-CoV-2 helicase formed up to three hydrogen bonds with 533 (Figure 12(e)).


### 2.7. Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA)

In this experiment, the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method was the utilized method to calculate the free binding energy of the SARS-CoV-2 helicase-533 complex. The MM-PBSA can evaluate the binding between a specific receptor and a ligand through the accurate calculation of the binding free energy of the ligand-protein complex. The MM-PBSA method utilizes both thermodynamic cycle and molecular dynamics (MD) methods to compute the binding free energy. The MM-PBSA calculates the binding free energies according to the following equation: ΔGbind. = Gcomp. − (Gprot. + Glig.).ΔGbind. refers to the total energy difference that was calculated as the difference between the energy at the bound-state (Gcomp.) and the sum of energy of both protein (Gprot.) and ligand (Glig.) before binding [98]. To compute the biding energy accurately, two main types of energies should be considered: first, the gas-phase interaction energy, which consists of van der Waals and electrostatic interactions; and, second, the solvation energy, which includes both polar and nonpolar components [99].The MM-PBSA, as a tool to calculate the free binding energies, has several advantages over other methods such as free energy perturbation and thermodynamic integration such as being faster, simpler, and producing consistent results with the experimental [100].The binding free energy of the SARS-CoV-2 helicase-533 complex was computed at the last stable 20 ns of the MD production run at a time interval of 100 ps from MD trajectories. The MM/PBSA method was utilized. Also, the MmPbSaStat.py script was employed to calculate the average free binding energy and its standard deviation/error from the output files that were obtained from g_mmpbsa. Compound 533, (Z)-6-(3-hydroxy-4-methoxystyryl)-4-methoxy-2H-pyran-2-one, showed a low binding free energy of −83 kJ/mol with the SARS-CoV-2 helicase (Figure 13(a)). The binding energy was stable during all the time of examination indicating the correct binding of the SARS-CoV-2 helicase-533 complex.Figure 13
MM-PBSA results of SARS-CoV-2 helicase-1552 complex.
(a)(b)


#### 2.7.1. Free Energy Decomposition

The total binding free energy of the SARS-CoV-2 helicase-533 complex was decomposed to analyze and understand the different components of that obtained binding energy as well as to disclose the contribution of each amino acid residue of the SARS-CoV-2 helicase in the binding with 533. The total binding free energy was decomposed into per amino acid residue contribution energy. This experiment gives a clearer idea about the essential amino acid residues that have favorable contributions to the binding process. It was found that THR-115, PHE-145, PRO-172, TYR-149, and LEU-411 residues of the protein contributed higher than −2 kJ/mol binding energy and thus they are crucial residues in the binding with 533 (Figure 13(b)).


## 2.1. Molecular Filtration Using Fingerprint Method

The ligand-basedin silico approach depends on the computation of chemical and physical properties of a molecule (ligand) and comparison of these properties with some biologically active compounds [48]. The fingerprint study is one of the ligand-based in silico methods that is vastly employed to predict the chemical structure’s similarity or dissimilarity of two compounds or more [49, 50]. During the fingerprint study, the computer converts the chemical descriptors of a molecule to mathematical symbols. The obtained data are presented as bit strings. These strings describe the presence (1) or absence (0) of a certain 2D fragment or atomic descriptor (property) in the test and reference molecules [51, 52]. The co-crystallized ligand is a molecule that has a very high affinity to bind to a specific protein forming a ligand-protein complex in a crystallized form [53]. In consequence, the chemical structure of the co-crystallized ligand could be utilized effectively as a starting point to design and discover a potential inhibitor against the target protein.Discovery studio 4.0 software was employed to examine the fingerprint similarity of 5956 natural compounds, which were derived from traditional Chinese medicine, againstVXG, the co-crystallized ligand of SARS-CoV-2 helicase (PDB ID: 5RMM). The experiment determined 187 compounds to be the most similar candidates to VXG (Table 1). The study compared the following descriptors (properties) in the chemical structures of the experiment set and VXG: H-bond acceptor [54], H-bond donor [55], charge [56], hybridization [57], positive ionizable atoms [58], negative ionizable atoms [59], halogens [60], aromatic groups [61], and aligning with the ALogP [62] of fragments as well as atoms.Table 1
Fingerprint similarity between 187 natural compounds andVXG.
Comp.SimilaritySASBSCComp.SimilaritySASBSCVXG11660054120.582205602580−1481381160100615850.581927586553−1321671169118−315660.581457439301152151180166−148100.581121197173−3125811531281354490.5808634312882334711431012323790.580247329113125364190−77651540.579775258−919637911431012332390.579154589563−13538011431012323560.57854523450−694111112425432580.57820545132634451186146−2035890.578089496404−424961196201−303620.577495272171825011163145354390.577398608599−1545071158136832550.577267452329252611431012315840.576874454333053311337033710.576792169127−35521173123−76230.576792169127−35541173121−7380.57608726561895551173121−754410.576077602591−1485771154901219820.575087494405−40610111639503980.57485961706191156131101140.5611111011465637116210848160.5611234547331169157−3920.556886931737561206193−401690.556886931737821220219−5451570.55670110828587921207178−414510.554878182162−167941196177−3051540.554307148101188031169147−37680.55418179157−138061210206−441900.55329910931578071183144−175570.55298167136−18181999676290.55176154−1023791213229−474950.54333316313434050.7771369304330.538462913754420.7381276397540.537764178165−12850.7311255415660.534743177165−112600.722117−4493540.534535178167−122080.69712412425970.53271171155−5910.68912211445960.53114816213941020.68912211446010.528302168152−21950.681109−6577520.527473192198−262800.6741166505680.526946176168−103440.65812118457900.526012182180−163700.658121184523680.523041227268−616720.65216790−14390.522222188194−22480.6491138531580.521978951671580.64911385351520.521898143108231000.646104−5622500.520408153128132460.644105−3615910.5197371581388840.638118194851530.517699117604951690.6361598471500.517241908768170.63312633402970.515789147119193960.62410646023700.515738213247−472450.622102−26451670.51567181185−15570.621065604130.51398614712019650.621065604910.5136991501261619420.614339497355−4351770.512195189203−2354150.614108592510−1384150.5119451501271654340.614017587502−13351510.511696175176−939410.6138334262402823670.511364225274−5940500.613551489343−355440.511299181188−1554410.613124626567−17251590.51111192147450390.613119458293−42660.510903164155240550.612984491347−373740.5105631451182135910.612751519393−654180.5104171471221923480.612622563465−10951550.510288124774239120.612319507374−532070.5102741491261735790.612128535420−8151490.51005203232−3739190.6114294282462651800.508197186200−2050550.611111440266142990.5070421441182251680.6115893851580.50632912071465390.609162100451760.505208194218−283420.60811015563970.505051150131165110.608178127−126140.505017151133157840.607193152−276220.504983152135145640.605182135−165120.50461516415922010.604186142−2051730.50358211253−455420.589905187151−216390.503356150132164920.58885169121−3500.50316515915072280.586319180141−141160.502982253337−8723650.584081521438−67800.502092120734632540.58408587551−13332880.502092240312−742060.583178624616−1704540.502075121754515710.582915464342−104190.5017421441212235910.582851503409−496210.501651521371415930.582512473358−192030.501449173179−7590.5823539946732910.501006249331−83640.58235399467860.515414212860.582278276201781180.51003466SA, bits number in both similar natural compounds and VXG; SB, bits number in the similar natural compounds but not in VXG; SC, bits number in VXG but not in the similar natural compounds.


## 2.2. Molecular Similarity

The difference between molecular similarity and fingerprint studies is that the fingerprint study computes the presence and/or absence of specific 2D atom paths and descriptors regarding fragments or substructures in the examined compounds [63]. Contrastingly, the molecular similarity calculates specific molecular descriptors considering the whole chemical structure of compounds. These descriptors are steric, topological, electronic, and/or physical [64]. Utilizing Discovery studio 4.0 software, a molecular similarity study was done on the most similar 187 compounds against VXG. The applied descriptors in this study (Figure 2 and Table 2) were partition coefficient (ALogp) [65], molecular weight (M. W) [66], H-bond donors (HBA) [67], H-bond acceptors (HBD) [68], rotatable bond numbers [69], number of rings as well as aromatic rings [70], and minimum distance [71] together with the molecular fractional polar surface area (MFPSA) [72]. The study was adapted to select the most similar 40 metabolites (Figure 3).Figure 2
Molecular similarity of the examined compounds andVXG.Table 2
Molecular descriptors of the examined 40 compounds andVXG.
Comp.ALogpM. WtHBAHBDRotatable bondsRingsAromatic ringsMFPSAMinimum distanceVXG0.71233.26312210.2374050.62219.24323210.3070.357842.12205.25313210.1870.4468162.1221.3322210.1810.4544331.69178.19313110.2390.5291001.09149.19212110.2520.5313422.76204.22312210.2240.5415331.97274.27514210.2290.5528181.53195.26323110.1870.5533641.44152.15312110.2850.5552461.91178.19313110.240.559590.8151.21222110.2640.56640.8151.21222110.2640.565391.57286.32423310.2390.5641501.54178.19322110.30.568852.03191.23201210.1430.5753981.81164.2213110.2020.5961951.45148.16202110.210.597571.23165.23223110.1650.614651.23165.23223110.1650.6148172.34205.3212210.1020.64851682.17276.37334210.1980.67251532.59249.35333210.1810.67851591.56195.26334110.2310.6911182.22203.28200210.1220.6953742.4317.36422310.2440.7062451.93162.19203110.1430.713481.77179.26213110.1050.715581.77179.26213110.1050.7151142.2194.23304110.1570.7316103.51230.26312210.1920.7383882.01333.36523310.2640.74351693.05260.37223210.1360.762601.93164.2204110.1410.765912.99188.22201210.1310.7691022.99188.22201210.1310.7692801.98177.24201210.0620.77551553.04263.38333210.1710.7826371.37336.77533310.2820.7842082.16191.27201210.0570.795800.9292.28631310.3370.81Figure 3
The most similar compounds toVXG.


## 2.3. Docking Studies

Computer-aided drug design applies various techniques to optimize natural products into potentially active against certain biological targets [44, 45]. A molecular docking technique was applied for 40 VXG similar compounds against SARS-CoV-2 helicase (PDB ID: 5RMM). The binding modes and affinities of these compounds were examined.The target protein was downloaded from Protein Data Bank (http://www.pdb.org), and Molecular Operating Environment (MOE .14) was used for the docking analysis. The docking process was validated by the re-docking of VXG inside the active pocket of the helicase protein. The root mean square deviation (RMSD) between the re-docked and the co-crystallized conformers was 0.74 A, which confirms the validity of the docking protocol (Figure 4).Figure 4
(a) 3D and (b) 2D superimposition of the re-docked conformer ofVXG over the co-crystallized one with an RMSD value of 0.74 A.
(a)(b)Out of the examined 40 compounds, 26 displayed correct binding modes as well as good free energies. The promising compounds are as follows: (+)-methylpseudoephedrine (48), (1R,2S)-ephedrine (57), (1R,2S)-N-methylephedrine (58), (1R,2S)-norephedrine (59), (3S)-2,2-dimethyl-3,5-dihydroxy-8-hydroxymethyl-3,4-dihydro-2H,6H-benzo[1,2-b5,4-b′]dipyran-6-one (80), 2-(4-(pyrrolidin-1-yl)phenyl)acetic acid (84), (4S,5R)-ephedroxane (85), 1-phenylpropane-1,2-dione (195), 2-methoxycinnamaldehyde (245), 2-methoxycinnamic acid (246), 2-methoxybenzoic acid (364), 3′-o-acetylhamaudol (374), 3′-o-propionylhamaudol (388), (R)-2-((R)-5-oxopyrrolidin-3-yl)-2-phenylacetic acid (405), (Z)-6-(3-hydroxy-4-methoxystyryl)-4-methoxy-2H-pyran-2-one (533), 6,7-dihydroxy-2-(2-phenylethyl)-5,6,7,8-tetrahydrochromone (539), 7-demethylsuberosin (610), 8-chloro-2-(2-phenylethyl)-5,6,7-trihydroxy-5,6,7,8-tetrahydrochromone (637), 3-((R)-hydroxy ((S)-1-methylpiperidin-2-yl)methyl)phenol (816), (R)-((S)-1-methylpiperidin-2-yl) (phenyl)methanol (817), 3-((1R,2S)-2-(dimethylamino)-1-hydroxypropyl)phenol (818), (R)-4-(1-Hydroxy-2-(methylamino)ethyl)-7,7-dimethyl-5,6,7,8-tetrahydronaphthalen-1-ol (5153),(R)-4-(1-hydroxy-2-(methylamino)ethyl)-8,8-dimethyl-6,7,8,9-tetrahydro-5H-benzo [7]annulen-1-ol (5155), (R)-2-ethyl-4-(1-hydroxy-2-(methylamino)ethyl)phenol (5159), (R)-2-((1S,2S,5S)-2-benzyl-5-hydroxy-4-methylcyclohex-3-en-1-yl)propane-1,2-diol (5168), and (1S,4R,5S)-4-benzyl-5-(2-hydroxypropan-2-yl)-2-methylcyclohex-2-en-1-ol (5169).The docking scores of the experienced ligands are predicted and summarized in Table3. The binding modes of the tested ligands inside the active site of the target protein were depicted. The binding poses of the top five compounds with the highest energy scores as well as the most perfect modes were selected for detailed discussion as representative examples.Table 3
∆G (in kcal/mole) of the most similar 40 compounds toVXG.
CompoundName∆G (kcal/mole)48(+)-Methylpseudoephedrine−15.9257(1R,2S)-Ephedrine−15.6458(1R,2S)-N-methylephedrine−15.2459(1R,2S)-Norephedrine−14.2964(1S,2S)-Norpseudoephedrine−13.6565(1S,2S)-Pseudoephedrine−11.6680(3S)-2,2-Dimethyl-3,5-dihydroxy-8-hydroxymethyl-3,4-dihydro-2H,6H-benzo[1,2-b5,4-b′]dipyran-6-one−14.27842-(4-(Pyrrolidin-1-yl) phenyl) acetic acid−15.4685(4S,5R)-Ephedroxane−14.0191(E)-3-Butylidene phthalide−13.31100(S)-Cathinone−12.6102(Z)-3-Butylidene phthalide−13.91141-(2,4-Dimethoxyphenyl)-1-propanone−13.531181,3,6,6-Tetramethyl-6,7-dihydroisoquinolin-8(5H)-one−12.41150(4S,5R)-Ephedroxane−13.021951-Phenylpropane-1,2-dione−11.072082,3,4-Trimethyl-5-phenyloxazolidine−13.762452-Methoxycinnamaldehyde−14.922462-Methoxycinnamic acid−14.222602-Phenylethyl acetate−13.292803,4-Dimethyl-5-phenyloxazolidine−13.213423-Butylidene-4-hydro-phthalide−13.243642-Methoxybenzoic acid−12.163743′-o-Acetylhamaudol−18.213883′-o-Propionylhamaudol−18.133984-(4-Hydroxyphenyl)-2-butanone−12.88405(R)-2-((R)-5-Oxopyrrolidin-3-yl)-2-phenylacetic acid−15.744334-Hydroxy-3-methoxycinnamaldehyde−13.11533(Z)-6-(3-Hydroxy-4-methoxystyryl)-4-methoxy-2H-pyran-2-one−17.15396,7-Dihydroxy-2-(2-phenylethyl)-5,6,7,8-tetrahydrochromone−18.366107-Demethylsuberosin−15.486378-Chloro-2-(2-phenylethyl)-5,6,7-trihydroxy-5,6,7,8-tetrahydrochromone−18.858163-((R)-Hydroxy ((S)-1-methylpiperidin-2-yl)methyl)phenol−17.09817(R)-((S)-1-Methylpiperidin-2-yl) (phenyl)methanol−15.958183-((1R,2S)-2-(Dimethylamino)-1-hydroxypropyl)phenol−15.795153(R)-4-(1-Hydroxy-2-(methylamino)ethyl)-7,7-dimethyl-5,6,7,8-tetrahydronaphthalen-1-ol−17.845155(R)-4-(1-Hydroxy-2-(methylamino)ethyl)-8,8-dimethyl-6,7,8,9-tetrahydro-5H-benzo[7]annulen-1-ol−16.595159(R)-2-Ethyl-4-(1-hydroxy-2-(methylamino)ethyl)phenol−17.895168(R)-2-((1S,2S,5S)-2-Benzyl-5-hydroxy-4-methylcyclohex-3-en-1-yl)propane-1,2-diol−18.45169(1S,4R,5S)-4-Benzyl-5-(2-hydroxypropan-2-yl)-2-methylcyclohex-2-en-1-ol−17.42VXG(3S,4R)-1-Acetyl-4-phenylpyrrolidine-3-carboxylic acid−19.37Starting with the binding interactions and orientation of the co-crystallized ligand (VXG) inside the active SARS-CoV-2 helicase (PDB ID: 5RMM), it revealed a binding affinity value of −19.37 kcal/mol. It showed a characteristic four hydrogen bonding interactions through carboxylate moiety of pyrrolidine ring with the essential amino acids SER486, ASN516, and ASN177. In addition, two hydrophobic interactions were formed between pyrrolidine ring and amino acid residues HIS554 and TYR515 (Figure 5).Figure 5
(a) 3D and (b) 2D ofVXG docked into the active site of SARS-CoV-2 helicase.
(a)(b)The results of docking studies showed that the tested ligands have orientations and binding interactions similar to that ofVXG against SARS-CoV-2 helicase. Three- and two-dimensional representations of binding modes of the most potent derivatives 533, 637, 84, 195, and 364 inside the active site of the target protein are depicted in Figures 6–10.Figure 6
(a) 3D and (b) 2D images of the docked compound533 into the active site of SARS-CoV-2 helicase.
(a)(b)Figure 7
(a) 3D and (b) 2D images of the docked compound637 into the active site of SARS-CoV-2 helicase.
(a)(b)Figure 8
(a) 3D and (b) 2D images of the docked compound84 into the active site of SARS-CoV-2 helicase.
(a)(b)Figure 9
(a) 3D and (b) 2D images of the docked compound195 into the active site of SARS-CoV-2 helicase.
(a)(b)Figure 10
(a) 3D and (b) 2D images of docked compound364 into the active site of SARS-CoV-2 helicase.
(a)(b)Compound533 exhibited an interesting binding mode similar to that of the co-crystallized ligand against SARS-CoV-2 helicase with a docking score of −17.10 kcal/mol. It keeps the hydrogen bonding interactions with the essential amino acids SER486, ASN516, and ASN177. Also, it formed two additional hydrogen bonds with ASN179 and SER485 residues. Furthermore, compound 533 was incorporated in two hydrophobic interactions with amino acid residues HIS554 and TYR515 (Figure 6).For compound637, the binding affinity was −18.85 kcal/mol. Such compound exhibited the best binding mode into the target protein, where it completely occupied the protein through seven hydrogen bonding interactions with SER486, ASN516, ASN177, SER485, and HIS554 residues. In addition, the terminal phenyl ring formed three hydrophobic interactions with HIS554, LEU412, and ARG560 residues (Figure 7).Concerning the binding mode of compound84 against SARS-CoV-2 helicase, the binding energy was −15.46 kcal/mol. The structural similarity of that compound with the co-crystallized ligand revealed the same binding mode against the receptor, where three hydrogen bonds were formed with SER486, ASN516, and ASN177 and three hydrophobic interactions were molded with HIS554 and PRO408 (Figure 8).The binding affinity of compound195 was −11.07 kcal/mol. Such affinity was represented by four hydrogen bonds with the key amino acids SER486, ASN516, and ASN177 and one hydrophobic interaction with HIS554 (Figure 9).Finally, analyzing the binding interactions of compound364 indicated a binding score of −12.16 kcal/mol. The carboxylate moiety formed three hydrogen bonding interactions with the key amino acids SER486, ASN516, and ASN177 while the phenyl ring was incorporated in hydrophobic interaction with HIS55 (Figure 10).


## 2.4. ADMET Studies

The likeness of any molecule to be approved as a drug depends greatly on its pharmacokinetic properties as well as its activity. Subsequently, the investigation of the ADMET profile of a molecule should be considered in the early stages of drug design and discovery to avoid the withdrawal possibility of the drug from the pharmaceutical market [73]. These descriptors identify the absorption, distribution, metabolism, excretion, as well as the toxicity of the examined compound. Although, there are different in vitro experiments that can determine the ADMET profile, in silico determination is an available and reliable tool with the profit of being faster, cheaper, as well as and lifesaver of the experimental animals [74].The predicted ADMET profiles of the 26 compounds that showed correct modes of binding besides Remdesivir, the reference drug, are shown in Table4 and Figure 11. Compounds 48, 58, 85, 245, 610, 816, 817, 5153, 5155, and 5169 were excluded because of their predicted strong ability to pass the blood-brain barriers which may be combined with a CNS toxicity. Among the excluded compounds, compounds 48, 58, 816, 817, 5153, and 5155 were predicted to be inhibitors against the CYP2D6 enzyme which would cause hepatotoxicity. The predicted intestinal absorption and aqueous solubility of all compounds were good to optimal.Table 4
ADMET profile of the 26 compounds with the best docking scores.
Comp.BBB levelaHIAbAqcCYP2D6dPPBe48104TF57204FF58104TF59304FF80303FF84203FT85103FF195203FT245103FT246203FT364204FT374203FF388303FF405304FF533303FT539303FT610103FT637303FT816203TF817103TF818204FF5153203TT5155203TT5159304FF5168204FT5169103FTSimeprevir432FTaBBB, ability to pass the blood-brain barrier, 1 is high, 2 is medium, 3 is low, and 4 is very low; bHIA, human intestinal absorption level, 0 is good, 1 is moderate, 2 is poor, and 3 is very poor; cAq, aqueous solubility level, 0 is extremely low, 1 is very low, 2 is low, 3 is good, and 4 is optimal; dCYP2D6, inhibition of CYP2D6 enzyme, T is an inhibitor and F is a noninhibitor; ePPB, F means less than 90% and T means more than 90%.Figure 11
Results of the ADMET study.


## 2.5. Toxicity Studies

The early prediction of toxicity is a crucial step that minimizes drug failure because of toxicity in the development stage or the clinical trials [75]. In silico prediction of toxicity is a credible approach that plays an essential role in drug design and discovery of lead compounds because in vitro and in vivo approaches are usually controlled by strict ethical regulations, time, and availability of resources [76], whereas the in silico prediction is based on a structure-activity relationship toxicology. The software compares the essential structural descriptors of the examined compounds with a huge library of hundreds of thousands of reported safe and toxic compounds [77] (Supporting data (available here)). Discovery studio 4.0 software was employed to predict the toxicity profile of the selected compounds after the ADMET study against 7 models. The applied models are FDA rat carcinogenicity [78, 79], mouse carcinogenic potency (TD50) [80], rat maximum tolerated dose (MTD) [81, 82], rat oral LD50 [83], rat chronic LOAEL [84, 85], ocular irritancy, and skin irritancy [86]. According to the obtained results (Table 5), compounds 80, 388, and 539 were eliminated due to the predicted high carcinogenic potency.Table 5
In silico toxicity profile of 15 compounds with good ADMET profile.
Comp.FDA rat carcinogenicity (mouse-female)TD50 (mg·kg−1·day−1)MTD (g·kg−1)Rat oral LD50 (g·kg−1)LOAEL (g·kg−1)Ocular irritancySkin irritancySimeprevirNot a carcinogen2.01380.0029670.2088350.0021057MildNone57Not a carcinogen306.6850.1269250.6785570.1091SevereNone59Not a carcinogen209.0130.1306691.178220.153373SevereNone80Multicarcinogen30.96540.3205630.1126450.0104398ModerateNone84Not a carcinogen86.27450.4320430.6515370.0466934SevereNone195Not a carcinogen734.3760.09202340.803940.58711MildNone246Not a carcinogen896.4370.1859080.9757830.0690491MildMild364Not a carcinogen1,152.330.1781251.100170.269177MildNone388Multicarcinogen77.44010.1582120.1285760.0111578SevereMild405Not a carcinogen1,019.870.3558430.5655660.0958513ModerateNone533Not a carcinogen587.5160.1090020.7653060.026065MildMild539Multicarcinogen71.15820.123540.3662150.0229232SevereNone637Not a carcinogen30.33650.1365950.4284340.0201137SevereMild818Not a carcinogen145.390.3557360.712450.0937066SevereNone5159Not a carcinogen227.5990.5599770.5267870.0808567SevereNone5168Not a carcinogen128.9110.2428542.096460.0494677ModerateMild


## 2.6. Molecular Dynamics (MD) Simulations

Despite the ability of molecular docking studies to expect the mode of binding of a compound inside a specific protein correctly, it has a serious defect that it deals with the protein as a rigid unit. Resultantly, it does not compute the conformational changes that happen in the protein because of ligand binding [87]. On the contrary, the MD simulations can adequately describe the behavior of a protein at the atomic level in full detail and at very accurate temporal resolution [88]. In accordance, the MD simulations have the advantage of being able to predict the conformational changes that occurred in the protein after ligand binding [89]. Furthermore, MD simulation studies can effectively compute various factors related to the energy of the protein-ligand complex for a determined time. Subsequently, it accurately describes the binding mode, stability, and flexibility of the ligand inside the target protein [90].The first successful MD simulation experiment of a protein (bovine pancreatic trypsin inhibitor) was published in Nature in 1977 [91]. Fortunately, because of the recently introduced supercomputer hardware, especially the advanced graphics processing units, MD simulation experiments became much more accessible, powerful, and accurate [92].In the MD simulation study, the forces on every atom of the examined ligand-protein complex are computed at every ultra-short time interval according to the basics of the “force field” [93]. The computed force field can be utilized to describe the position and velocity of atoms at each time interval. The force field is a physical expression that describes the functional potential energy of atoms. The force field is calculated based on Newton’s laws of motion considering bonded interactions (bonds, angles, and dihedrals) in addition to nonbonded interactions (van der Waals potentials and Coulomb potentials) between all atoms of the complex. This step is repeated billions of times to produce the atomic trajectories for a specific time interval [94].Several MD simulation studies were employed to investigate the stability and mimic the dynamic of compound533, (Z)-6-(3-hydroxy-4-methoxystyryl)-4-methoxy-2H-pyran-2-one, that exhibited the best docking score inside SARS-CoV-2 helicase for 100 ns.First, the interaction of a ligand inside the active site of a protein leads to some changes in the structure of that protein [95]. In consequence, the dynamics and the conformational changes of the SARS-CoV-2 helicase-533 complex were computed as root mean square deviation (RMSD) to detect the stability due to binding. It is observed that the SARS-CoV-2 helicase and 533 exhibited lower RMSD with no major fluctuations indicating their greater stability (Figure 12(a)). Interestingly, the SARS-CoV-2 helicase-533 complex was stable till 90 ns∼. Although the SARS-CoV-2 helicase-533 complex showed a minor fluctuation later, it reached equilibrium again.Figure 12
M D simulation results: (a) RMSD values of SARS-CoV-2 helicase,533, and SARS-CoV-2 helicase-533 complex; (b) RMSF for SARS-CoV-2 helicase; (c) Rg of SARS-CoV-2 helicase; (d) SASA of SARS-CoV-2 helicase; and (e) H-bonding between SARS-CoV-2 helicase-533 complex in the MD run.
(a)(b)(c)(d)(e)Second, the flexibility of SARS-CoV-2 helicase was calculated in terms of root mean square fluctuation (RMSF) to calculate the differences in flexibility in the SARS-CoV-2 helicase-533 complex during the 100 ns of the MD simulations. The decrease of RMSF values during the MD simulation in 50–100 residue areas (Figure 12(b)) denotes that SARS-CoV-2 residues were more rigid and stabilized after binding to 533.Third, the radius of gyration (Rg) is a crucial parameter that is linked to the protein stability according to the change in its volume. Rg is defined as the root mean square distance (RMSD) of a weighted mass group of atoms from their mass center [96, 97]. Thus, the calculation of Rg identifies the dimensions as well as the compactness of the SARS-CoV-2 helicase-533 complex. The lower degree of fluctuation throughout the simulation period indicates the greater compactness of a system. The Rg of the SARS-CoV-2 helicase-533 complex was found to be lower than the starting period (Figure 12(c)) displaying compactness and stability.Fourth, the interaction between protein-ligand complexes and solvents was measured by solvent accessible surface area (SASA) over the simulation period. So, the SASA of the SARS-CoV-2 helicase-533 complex was calculated to provide the extent of the conformational changes that occurred during binding. Interestingly, SARS-CoV-2 helicase featured a reduction of the surface area showing a relatively lower SASA value than the starting period (Figure 12(d)).Finally, hydrogen bonding between a SARS-CoV-2 helicase-533 complex is essential to stabilize the structure. MD simulation studies showed that the highest number of conformations of the SARS-CoV-2 helicase formed up to three hydrogen bonds with 533 (Figure 12(e)).


## 2.7. Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA)

In this experiment, the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method was the utilized method to calculate the free binding energy of the SARS-CoV-2 helicase-533 complex. The MM-PBSA can evaluate the binding between a specific receptor and a ligand through the accurate calculation of the binding free energy of the ligand-protein complex. The MM-PBSA method utilizes both thermodynamic cycle and molecular dynamics (MD) methods to compute the binding free energy. The MM-PBSA calculates the binding free energies according to the following equation: ΔGbind. = Gcomp. − (Gprot. + Glig.).ΔGbind. refers to the total energy difference that was calculated as the difference between the energy at the bound-state (Gcomp.) and the sum of energy of both protein (Gprot.) and ligand (Glig.) before binding [98]. To compute the biding energy accurately, two main types of energies should be considered: first, the gas-phase interaction energy, which consists of van der Waals and electrostatic interactions; and, second, the solvation energy, which includes both polar and nonpolar components [99].The MM-PBSA, as a tool to calculate the free binding energies, has several advantages over other methods such as free energy perturbation and thermodynamic integration such as being faster, simpler, and producing consistent results with the experimental [100].The binding free energy of the SARS-CoV-2 helicase-533 complex was computed at the last stable 20 ns of the MD production run at a time interval of 100 ps from MD trajectories. The MM/PBSA method was utilized. Also, the MmPbSaStat.py script was employed to calculate the average free binding energy and its standard deviation/error from the output files that were obtained from g_mmpbsa. Compound 533, (Z)-6-(3-hydroxy-4-methoxystyryl)-4-methoxy-2H-pyran-2-one, showed a low binding free energy of −83 kJ/mol with the SARS-CoV-2 helicase (Figure 13(a)). The binding energy was stable during all the time of examination indicating the correct binding of the SARS-CoV-2 helicase-533 complex.Figure 13
MM-PBSA results of SARS-CoV-2 helicase-1552 complex.
(a)(b)


### 2.7.1. Free Energy Decomposition

The total binding free energy of the SARS-CoV-2 helicase-533 complex was decomposed to analyze and understand the different components of that obtained binding energy as well as to disclose the contribution of each amino acid residue of the SARS-CoV-2 helicase in the binding with 533. The total binding free energy was decomposed into per amino acid residue contribution energy. This experiment gives a clearer idea about the essential amino acid residues that have favorable contributions to the binding process. It was found that THR-115, PHE-145, PRO-172, TYR-149, and LEU-411 residues of the protein contributed higher than −2 kJ/mol binding energy and thus they are crucial residues in the binding with 533 (Figure 13(b)).


## 2.7.1. Free Energy Decomposition

The total binding free energy of the SARS-CoV-2 helicase-533 complex was decomposed to analyze and understand the different components of that obtained binding energy as well as to disclose the contribution of each amino acid residue of the SARS-CoV-2 helicase in the binding with 533. The total binding free energy was decomposed into per amino acid residue contribution energy. This experiment gives a clearer idea about the essential amino acid residues that have favorable contributions to the binding process. It was found that THR-115, PHE-145, PRO-172, TYR-149, and LEU-411 residues of the protein contributed higher than −2 kJ/mol binding energy and thus they are crucial residues in the binding with 533 (Figure 13(b)).


## 3. Methods

### 3.1. Molecular Similarity

Discovery studio 4.0 software was used [101, 102] (see Section 3 in Supplementary data).


### 3.2. Fingerprints Studies

Discovery studio 4.0 software was used [103, 104] (see Section 3 in Supplementary data).


### 3.3. Docking Studies

Docking studies were done against the target enzyme using Discovery studio 4.0 software [105, 106] (see Section 3 in Supplementary data).


### 3.4. ADMET Analysis

Discovery studio 4.0 was used [40, 107] (see Section 3 in Supplementary data).


### 3.5. Toxicity Studies

Discovery studio 4.0 software was used [108–110] (see Section 3 in Supplementary data).


### 3.6. Molecular Dynamics Simulation

The system was prepared using the web-based CHARMM-GUI [111–113] utilizing CHARMM36 force field [114] and NAMD 2.13 [115] package. The TIP3P explicit solvation model was used (see supporting data (available here)).


### 3.7. MM-PBSA Studies

The g_mmpbsa package of GROMACS was utilized to calculate the MM/PBSA (see supporting data (availablehere)).


## 3.1. Molecular Similarity

Discovery studio 4.0 software was used [101, 102] (see Section 3 in Supplementary data).


## 3.2. Fingerprints Studies

Discovery studio 4.0 software was used [103, 104] (see Section 3 in Supplementary data).


## 3.3. Docking Studies

Docking studies were done against the target enzyme using Discovery studio 4.0 software [105, 106] (see Section 3 in Supplementary data).


## 3.4. ADMET Analysis

Discovery studio 4.0 was used [40, 107] (see Section 3 in Supplementary data).


## 3.5. Toxicity Studies

Discovery studio 4.0 software was used [108–110] (see Section 3 in Supplementary data).


## 3.6. Molecular Dynamics Simulation

The system was prepared using the web-based CHARMM-GUI [111–113] utilizing CHARMM36 force field [114] and NAMD 2.13 [115] package. The TIP3P explicit solvation model was used (see supporting data (available here)).


## 3.7. MM-PBSA Studies

The g_mmpbsa package of GROMACS was utilized to calculate the MM/PBSA (see supporting data (availablehere)).


## 4. Conclusion

Twelve of 5956 TCM compounds were suggested to be the potential inhibitors against SARS-CoV-2 helicase (PDB ID: 5RMM). The compounds were selected according to structural similarity and fingerprint studies withVXG, the co-crystallized ligand of the target protein. Then, molecular docking studies were carried out. Then, ADMET and toxicity studies were preceded to select the following metabolites: (1R,2S)-ephedrine (57), (1R,2S)-norephedrine (59), 2-(4-(pyrrolidin-1-yl)phenyl)acetic acid (84), 1-phenylpropane-1,2-dione (195), 2-methoxycinnamic acid (246), 2-methoxybenzoic acid (364), (R)-2-((R)-5-oxopyrrolidin-3-yl)-2-phenylacetic acid (405), (Z)-6-(3-hydroxy-4-methoxystyryl)-4-methoxy-2H-pyran-2-one (533), 8-chloro-2-(2-phenylethyl)-5,6,7-trihydroxy-5,6,7,8-tetrahydrochromone (637), 3-((1R,2S)-2-(dimethylamino)-1-hydroxypropyl)phenol (818), (R)-2-ethyl-4-(1-hydroxy-2-(methylamino)ethyl)phenol (5159), and (R)-2-((1S,2S,5S)-2-benzyl-5-hydroxy-4-methylcyclohex-3-en-1-yl)propane-1,2-diol (5168). Among them, compounds 84, 195, 364, 533, and 637 showed the best docking scores. Interestingly, compound 533, the one with the highest docking score, bonded favorably to the target protein with low energy and optimum dynamics according to advanced MD simulation studies over 100 ns.



---

*Source: 7270094-2022-08-22.xml*