key: cord-1029415-e7cgfhe4
authors: Ibrahim, Mahmoud A. A.; Abdelrahman, Alaa H. M.; Mohamed, Tarik A.; Atia, Mohamed A. M.; Al-Hammady, Montaser A. M.; Abdeljawaad, Khlood A. A.; Elkady, Eman M.; Moustafa, Mahmoud F.; Alrumaihi, Faris; Allemailem, Khaled S.; El-Seedi, Hesham R.; Paré, Paul W.; Efferth, Thomas; Hegazy, Mohamed-Elamir F.
title: In Silico Mining of Terpenes from Red-Sea Invertebrates for SARS-CoV-2 Main Protease (M(pro)) Inhibitors
date: 2021-04-05
journal: Molecules
DOI: 10.3390/molecules26072082
sha: 3db8a1b5d80f30d408cffe17d45482d6d044ba57
doc_id: 1029415
cord_uid: e7cgfhe4

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the causative agent for the COVID-19 pandemic, which generated more than 1.82 million deaths in 2020 alone, in addition to 83.8 million infections. Currently, there is no antiviral medication to treat COVID-19. In the search for drug leads, marine-derived metabolites are reported here as prospective SARS-CoV-2 inhibitors. Two hundred and twenty-seven terpene natural products isolated from the biodiverse Red-Sea ecosystem were screened for inhibitor activity against the SARS-CoV-2 main protease (M(pro)) using molecular docking and molecular dynamics (MD) simulations combined with molecular mechanics/generalized Born surface area binding energy calculations. On the basis of in silico analyses, six terpenes demonstrated high potency as M(pro) inhibitors with ΔG(binding) ≤ −40.0 kcal/mol. The stability and binding affinity of the most potent metabolite, erylosides B, were compared to the human immunodeficiency virus protease inhibitor, lopinavir. Erylosides B showed greater binding affinity towards SARS-CoV-2 M(pro) than lopinavir over 100 ns with ΔG(binding) values of −51.9 vs. −33.6 kcal/mol, respectively. Protein–protein interactions indicate that erylosides B biochemical signaling shares gene components that mediate severe acute respiratory syndrome diseases, including the cytokine- and immune-signaling components BCL2L1, IL2, and PRKC. Pathway enrichment analysis and Boolean network modeling were performed towards a deep dissection and mining of the erylosides B target–function interactions. The current study identifies erylosides B as a promising anti-COVID-19 drug lead that warrants further in vitro and in vivo testing.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) belongs to the family Coronaviridae (genus Betacoronavirus; subgenus Sarbecovirus) with a viral envelope encapsulating a single-stranded RNA genome. While most coronaviruses infect mammals and birds, the animal host responsible for the transmission of SARS-CoV-2 to men has yet to be identified [1] [2] [3] .

The SARS-CoV-2 genome contains 11 genes encoding 29 proteins and peptides ( www.ncbi.nlm.nih.gov/nuccore/NC_045512.2?report=graph). Four proteins constitute the viral structure, including the spike or S protein responsible for binding with the host cell receptor [4] . The S protein binds to the angiotensin-converting enzyme 2 (ACE2), which is a necessary step for viral entry into the cell. Two currently employed anti-SARS-CoV-2 vaccines utilize nanoparticle-encapsulated mRNA coding for the spike protein.

These vaccines activate human immune responses to target the S protein, and about 95% of vaccinated people are protected from severe courses of the Coronavirus disease 2019 (COVID-19) [5] . In the case of more contagious variants such as a strain originating from Brazil (known as P.1), the genome has 17 unique mutations, including three in the receptor-binding domain of the spike protein. Mutations in the S protein are thought to be responsible for the variable ability of the spike protein to bind to ACE2, and specifically, the P.1 mutant reveals increased affinity and penetration into the host. Recently, an increasing number of escape mutations in the spike protein threaten the effectiveness of vaccines [6] .

In seeking alternative approaches to vaccination, other SARS-CoV-2 proteins can be investigated that control the virus replication. Essential nonstructural proteins are initially expressed as two large polyproteins that are then processed into 16 peptide components. The main protease (M pro ) initially cleaves the polyproteins into 11 fragments. Inhibitors that block the catalytic M pro activity effectively interrupt the viral replication. One such inhibitor has been recently identified [7] . The crystal structure of SARS-CoV-2 main protease provides a basis for the design of improved alpha-ketoamide inhibitors [7] and M pro appears to be a promising target for designing novel small molecule inhibitors. Since that original study, several promising inhibitor leads have been identified using in silico enzyme modeling [8] . Structure-based computational modeling of ligand-receptor interactions was used by Ibrahim et al. to identify potential M pro inhibitors [9] [10] [11] [12] [13] . Natural products hold a vital role in discovering novel and effective therapeutics to combat the present COVID-19 pandemic. Among natural products, flavonoids, alkaloids, and terpenoids have attracted great attention as prospective SARS-CoV-2 inhibitors [14] [15] [16] . Recognizing that marine invertebrates are promising organisms for biologically active metabolites including anti-inflammatory, antibacterial, antifungal, antimalarial, antitumor, and antiviral activity [17, 18] , here biologically active terpene metabolites identified from a coral reef community unique to the Red Sea [19] were screened for in silico binding affinities against SARS-CoV-2 M pro . Previously characterized metabolites from this natural-product pool include alismol and aromadendrane sesquiterpenes derived from Litophyton arboretum [20] that exhibit inhibitory activity against the HIV-1 protease (HIV-1 PR) (IC 50 7 µM); palustrol, a sesquiterpene from Sarcophyton trocheliophorum that has antibacterial activity (MIC 6.6-11.1 µM) [21] ; and 12(S)-Hydroperoxylsarcoph-10-ene, a cembrane diterpene from Sarcophyton glaucum that was reported to exhibit potent anticancer activity via the inhibition of Cyp1A activity (p < 0.01) with IC 50 values of 2.7 nM [22] . On the basis of the predicted docking scores, the most potent inhibitors are submitted to molecular dynamics (MD) simulations combined with binding energy calculations using a molecular mechanics/generalized Born surface area approach.

Since the main protease (M pro ) of SARS-CoV-2 plays an indispensable role in viral reproduction, small molecules were screened based on in silico molecular docking calculations and MD simulations for prospective M pro inhibitors. Marine natural products identified from the Red Sea provided the source for metabolite screening.

Two hundred and twenty-seven terpene natural products isolated from the biodiverse Red-Sea ecosystem were screened against the SARS-CoV-2 main protease (M pro ) using molecular docking technique. Molecular docking calculations resulted in 27 of the screened compounds exhibiting a higher binding affinity than lopinavir: an inhibitor of SARS-CoV-2 main protease (M pro ) that was proposed as a treatment for COVID-19 on the basis of in vitro activity, preclinical studies, and observational studies [23] . While docking scores ranged from −4.3 to −12.3 kcal/mol, 12% of the compounds scored below −9.8 kcal/mol (Table S1 ). AutoDock4.2.6 software was utilized to carry out all molecular docking calculations. Binding affinities, 2D chemical structures, and features of the 27 most promising natural products towards SARS-CoV-2 M pro are summarized in Table 1 . 2D docking positions with proximal amino acid residues within the M pro active site are depicted in Figure S1 . Most of these compounds demonstrate similar M pro binding modes within the binding pocket, forming hydrogen bonds with CYS145, HIS164, and GLU166, which can account for the high binding affinities (Table 1 and Figure S1 ). The 2D and 3D representations of the interactions of the top three potent marine natural products (MNPs) and lopinavir with key amino acid residues of SARS-CoV-2 M pro are depicted in Figure 1 and Figure S2 , respectively.

Depresosterol (190) isolated from Lobophytum depressum, exhibited the highest binding affinity of the compounds screened against SARS-CoV-2 M pro , with a docking score of −12.3 kcal/mol. in silico M pro binding in the active site indicated that the methanolic hydroxyl group exhibited two hydrogen bonds with a backbone carboxylate of GLU166 with bond lengths of 1.99 and 2.55 Å, respectively ( Figure 1 and Table 1 ). In addition, the hydroxyl unit of 2-methylpropan-2-ol affords three hydrogen bonds with a backbone NH and carbonyl group of ASN142 with bond lengths of 2.24, 2.68, and 2.04 Å, respectively ( Figure 1 and Table 1 ). Moreover, the hydroxy group of 2-propanol exhibited a hydrogen bond with the backbone carbonyl group of ASN142 with a bond length of 1.96 Å (Figure 1 , Figure S2 and Table 1 ). The oxygen of the oxirane ring interacted with the backbone imidazole ring of HIS41, and the thiol group of CYS145 with bond lengths of 2.17 and 2.70 Å, respectively ( Figure 1 and Table 1 ). The hydroxy group of the cyclohexanol ring contributed two hydrogen bonds with NH and the carbonyl group of TYR26 with bond lengths of 2.15 and 2.66 Å, respectively ( Figure 1 and Table 1) .

3β-25-Dihydroxy-4-methyl-5α,8α-epidioxy-2-ketoergost-9-ene (178) isolated from Sinularia candidula, exhibited the second highest binding affinity of the compounds screened against SARS-CoV-2 M pro with a docking score of −12.2 kcal/mol. in silico M pro binding in the active site indicated that the hydroxy group of the hydroxycyclohexanone ring participates in four hydrogen bonds with the backbone carbonyl of LEU141, OH and NH of SER144, and NH of CYS145 with bond lengths of 2.08, 1.97, 2.28, and 2.49 Å, respectively ( Figure 1 and Table 1 ). Moreover, the carbonyl group of the hydroxycyclohexanone ring demonstrates two hydrogen bonds with an imidazole ring of HIS163 and backbone OH of SER144 with bond lengths of 2.08 and 2.66 Å, respectively ( Figure 1 and Table 1 ). The hydroxy group of 2-methylpropan-2-ol displays two hydrogen bonds with a backbone carbonyl group of THR190 and NH of GLN192 with bond lengths of 1.80 and 2.26 Å, respectively ( Figure 1 and Table 1 ). The oxygen atom of the methoxy group is also hydrogen bound to the NH of ASN142 with a bond length of 2.96 Å ( Figure 1 and Table 1 ). Compared with 190, 178, and 226, lopinavir presented a similar binding affinity towards M pro with a docking score of −9.8 kcal/mol. Scrutinizing the binding mode of lopinavir inside the active site of M pro revealed that the NH of the tetrahydro-1-methylpyrimidin-2(1H)-one ring exhibited two hydrogen bonds with a backbone hydroxy group of SER144 and a carbonyl group of LEU141 with bond lengths of 3.09 and 1.96 Å, respectively. lopinavir presented a similar binding affinity towards M pro with a docking score of −9.8 kcal/mol. Scrutinizing the binding mode of lopinavir inside the active site of M pro revealed that the NH of the tetrahydro-1-methylpyrimidin-2(1H)-one ring exhibited two hydrogen bonds with a backbone hydroxy group of SER144 and a carbonyl group of LEU141 with bond lengths of 3.09 and 1.96 Å, respectively. kcal/mol. Scrutinizing the binding mode of lopinavir inside the active site of M pro revealed that the NH of the tetrahydro-1-methylpyrimidin-2(1H)-one ring exhibited two hydrogen bonds with a backbone hydroxy group of SER144 and a carbonyl group of LEU141 with bond lengths of 3.09 and 1.96 Å, respectively. Erylosides B (226) isolated from Erylus lendenfeldi also exhibited a high binding affinity against SARS-CoV-2 M pro with a docking score of −12.1 kcal/mol. in silico M pro binding in the active site indicated that the five carbonyl groups of methyl acetates exhibit seven hydrogen bonds with a backbone OH of TYR54, imidazole ring of HIS41, SH of CYS145, imidazole ring of HIS163, NH of ASN142, NH of GLY143, and NH 2 of GLN189 with bond lengths of 3.05, 1.96, 2.04, 1.93, 1.75, 2.82, and 1.91 Å, respectively ( Figure 1 and Table 1 ). Moreover, the oxygen of the (cyclohexyloxy) cyclohexane ring forms a hydrogen bond with the backbone NH of GLU166 with a bond length of 2.37 Å (Figure 1 and Table 1 ). The oxygen atom of the methoxy group is also hydrogen bound to the NH of ASN142 with a bond length of 2.96 Å (Figure 1 and Table 1 ). Compared with 190, 178, and 226, lopinavir presented a similar binding affinity towards M pro with a docking score of −9.8 kcal/mol. Scrutinizing the binding mode of lopinavir inside the active site of M pro revealed that the NH of the tetrahydro-1-methylpyrimidin-2(1H)-one ring exhibited two hydrogen bonds with a backbone hydroxy group of SER144 and a carbonyl group of LEU141 with bond lengths of 3.09 and 1.96 Å, respectively.

The carbonyl group of the tetrahydro-1-methylpyrimidin-2(1H)-one ring was also observed to form a hydrogen bond with the backbone NH of GLY143 with a bond length of 2.01 Å, and the NH of the methylacetamide group participates in a hydrogen bond with the carbonyl group of HIS164 with a bond length of 2.62 Å (Figure 1 and Table 1 ). The current results provide quantitative data of the binding affinities of 190, 178, 226, and lopinavir as promising SARS-CoV-2 M pro inhibitors.

MD simulations probe the stability of the ligand-enzyme complexes, conformational flexibilities, structural details, and the dependability of ligand-enzyme affinities [24, 25] . Therefore, those natural products with low M pro docking scores were submitted for MD simulations followed by binding energy calculations. The simulations were conducted with an implicit water solvent for 250 ps, and the molecular mechanics/generalized Born surface area (MM/GBSA) approach was applied to estimate the corresponding binding affinities; thereby diminishing the time and computational costs (see computational methodology section for details). The calculated MM/GBSA binding energies for the pre-screened natural products are summarized in Table S2 . Six compounds had lower binding energies (∆G binding ) than lopinavir (calc. −39.4 kcal/mol). These metabolites were further subjected to a 10 ns MD simulation in an explicit water solvent to obtain more accurate M pro binding affinities. MD/MM/GBSA binding energies were estimated (Figure 2 ). Four of these compounds exhibited lower binding energies (∆G binding ) compared to lopinavir (calc. No appreciable variance between the calculated MM/GBSA binding energy for the 226-M pro complex throughout the 50 ns MD simulation and the corresponding MM/GBSA binding energy throughout the 100 ns MD simulation were observed (Figure 2) . Compared with lopinavir, 226 showed a higher binding affinity towards M pro over the 100 ns MD simulation with an average ∆G binding of −51.9 kcal/mol. Hydrogen bonding, pi-based, hydrophobic, and van der Waals interactions with essential amino acid residues within the M pro active site are all thought to contribute to the strong binding constant. The 2D and 3D representations for 226 and lopinavir within the M pro active site over 100 ns are shown in Figure 3 and S3, respectively. Interestingly, 226 maintained hydrogen bonding with the proximal amino acid residues of M pro over the 100 ns MD course (Figure 3 ). Structural insights into the binding mode of the 226 with the M pro demonstrated that the carbonyl groups of methyl acetates form four hydrogen bonds with the backbone imidazole of the HIS41, OH of SER146, the carbonyl group of ASN142, and NH2 of GLN189, with bond lengths of 2.72, 2.76, 3.37, and 1.89 Å, respectively ( Figure 3) ; while the oxygen of the (cyclohexyloxy) cyclohexane ring forms hydrogen bonds with the backbone NH2 of GLN189 with a bond length of 2.62 Å (Figure 3) . Furthermore, the oxygen of tetrahydro-2Hpyran bonds with the backbone NH2 of GLN189 with a bond length of 2.06 Å (Figure 3 ). of these compounds exhibited lower binding energies (ΔGbinding) compared to lopinavir (calc. −35.4 kcal/mol). Moreover, those potent MNPs were chosen and submitted for 50 ns MD simulations in the explicit water solvent, and the corresponding binding energies were evaluated (Figure 2 ). Only erylosides B (226) exhibited steady binding throughout the simulation, while 202, 151, and 224 showed an increase in MM/GBSA binding energies over the simulated time. For instance, the calculated MM/GBSA binding energies for 226 against M pro were −50.8, −48.8 and −50.0 kcal/mol over 250 ps implicit-solvent MD, 10 ns explicit-solvent MD, and 50 ns explicit-solvent MD simulations, respectively. This shows the importance of long MD simulations to predict metabolite-M pro binding affinities. Therefore, MD simulations for the 226-M pro complex were prolonged to 100 ns, and the corresponding MM/GBSA binding energy was calculated (Figure 2 ). No appreciable variance between the calculated MM/GBSA binding energy for the 226-M pro complex throughout the 50 ns MD simulation and the corresponding MM/GBSA binding energy throughout the 100 ns MD simulation were observed (Figure 2) . Compared with lopinavir, 226 showed a higher binding affinity towards M pro over the 100 ns MD simulation with an average ΔGbinding of −51.9 kcal/mol. Hydrogen bonding, pibased, hydrophobic, and van der Waals interactions with essential amino acid residues within the M pro active site are all thought to contribute to the strong binding constant. The 2D and 3D representations for 226 and lopinavir within the M pro active site over 100 ns are shown in Figure 3 and S3, respectively. Interestingly, 226 maintained hydrogen bonding with the proximal amino acid residues of M pro over the 100 ns MD course (Figure 3 ). Structural insights into the binding mode of the 226 with the M pro demonstrated that the carbonyl groups of methyl acetates form four hydrogen bonds with the backbone imidazole of the HIS41, OH of SER146, the carbonyl group of ASN142, and NH2 of GLN189, with bond lengths of 2.72, 2.76, 3.37, and 1.89 Å, respectively ( Figure 3) ; while the oxygen of the (cyclohexyloxy) cyclohexane ring forms hydrogen bonds with the backbone NH2 of GLN189 with a bond length of 2.62 Å ( Figure 3) . Furthermore, the oxygen of tetrahydro-2H-pyran bonds with the backbone NH2 of GLN189 with a bond length of 2.06 Å (Figure 3 ).

While lopinavir also displayed high M pro binding energy over the 100 ns MD simulation at ΔGbinding −33.6 kcal/mol, only three hydrogen bonds were observed with the proximal amino acid residues of M pro (Figure 3) . A comparison of the erylosides B (226) and lopinavir revealed that the binding affinity of erylosides B (226) was approximately two times higher than that of lopinavir. While lopinavir also displayed high M pro binding energy over the 100 ns MD simulation at ∆G binding −33.6 kcal/mol, only three hydrogen bonds were observed with the proximal amino acid residues of M pro (Figure 3) . A comparison of the erylosides B (226) and lopinavir revealed that the binding affinity of erylosides B (226) was approximately two times higher than that of lopinavir.

The estimated MM/GBSA binding energies were further decomposed into separate components to indicate the force in the binding of M pro with erylosides B (226) and lopinavir ( Table 2 ). E vdw was a significant contributor to the erylosides B (226)-and lopinavir-M pro binding affinities with average values of −71.2 and −45.6 kcal/mol, respectively. E ele was efficient with an average value of −30.5 and −22.1 kcal/mol for the erylosides B (226)-and lopinavir-M pro binding affinities, respectively. The binding energies of erylosides B (226)-and lopinavir-M pro complexes were further decomposed at the per-residue level and the residues with absolute energy contribution < −0.50 kcal/mol were shown ( Figure 4) . MET165, GLU166, PRO168, and GLN189 in the M pro complex favorably participate with erylosides B (226) and lopinavir. There was a considerable contribution by GLN189 to the total binding free energy with values of −4.6 and −3.3 kcal/mol for erylosides B (226) and lopinavir, respectively. Additionally, the hydrophobic residues share greater energies as a result of the formation of the hydrophobic interactions between erylosides B (226) and hydrophobic residues. The estimated MM/GBSA binding energies were further decomposed into separate components to indicate the force in the binding of M pro with erylosides B (226) and lopinavir ( Table 2) The binding energies of erylosides B (226)-and lopinavir-M pro complexes were further decomposed at the per-residue level and the residues with absolute energy contribution < −0.50 kcal/mol were shown ( Figure 4) . MET165, GLU166, PRO168, and GLN189 in the M pro complex favorably participate with erylosides B (226) and lopinavir. There was a considerable contribution by GLN189 to the total binding free energy with values of −4.6 and −3.3 kcal/mol for erylosides B (226) and lopinavir, respectively. Additionally, the hydrophobic residues share greater energies as a result of the formation of the hydrophobic interactions between erylosides B (226) and hydrophobic residues. 

To further confirmed the stability and behavior of erylosides B (226) complexed with M pro , structural and energetic analyses were executed during 100 ns MD simulations and compared to those of lopinavir. Control of the structural stability of the studied system was accomplished by investigating root-mean-square deviation (RMSD), center-of-mass (CoM) distance, hydrogen bond length, and binding energy per frame. 

To further confirmed the stability and behavior of erylosides B (226) complexed with M pro , structural and energetic analyses were executed during 100 ns MD simulations and compared to those of lopinavir. Control of the structural stability of the studied system was accomplished by investigating root-mean-square deviation (RMSD), center-of-mass (CoM) distance, hydrogen bond length, and binding energy per frame.

The global structural stability of erylosides B (226) and lopinavir in complex with M pro was estimated over the 100 ns MD simulations by measuring the correlation between the binding energy per frame and time ( Figure 5 ). Overall stabilities for erylosides B (226) and lopinavir exhibited average binding energies (∆G binding ) of −51.9 and −33.6 kcal/mol, respectively. According to this analysis, all complexes conserved stability throughout the 100 ns MD simulation. The global structural stability of erylosides B (226) and lopinavir in complex with M pro was estimated over the 100 ns MD simulations by measuring the correlation between the binding energy per frame and time ( Figure 5 ). Overall stabilities for erylosides B (226) and lopinavir exhibited average binding energies (ΔGbinding) of −51.9 and −33.6 kcal/mol, respectively. According to this analysis, all complexes conserved stability throughout the 100 ns MD simulation. 

Hydrogen bond analyses for erylosides B (226)-and lopinavir-M pro complexes over a 100 ns MD simulation were performed (Table 3) . Four hydrogen bonds with GLN189, GLU166, CYS145, and ASN142 were observed within the erylosides B (226)-M pro complex. Average bond lengths were 2.9, 2.8, 2.9, and 2.7 Å, with occupation percentages of 96, 92, 91, and 83%, respectively. In contrast, two hydrogen bonds were detected for the lopinavir-M pro complex, with GLN189 and GLY143 showing an average bond length of 2.8 and 2.7 Å and an occupation percentage of 85.6 and 75.6%, respectively. Overall, these hydrogen bond analyses supported the finding of high stability for the erylosides B (226)-M pro complex compared to lopinavir-M pro . Lopinavir @O12-H28 2.7 158 75.6 a The hydrogen bonds are inspected by the acceptor-H-donor angle of >120° and acceptor-donor atom distance < 3.5 Å. b Occupancy is employed to estimate the strength and stability of the hydrogen bond. 

Hydrogen bond analyses for erylosides B (226)-and lopinavir-M pro complexes over a 100 ns MD simulation were performed (Table 3) . Four hydrogen bonds with GLN189, GLU166, CYS145, and ASN142 were observed within the erylosides B (226)-M pro complex. Average bond lengths were 2.9, 2.8, 2.9, and 2.7 Å, with occupation percentages of 96, 92, 91, and 83%, respectively. In contrast, two hydrogen bonds were detected for the lopinavir-M pro complex, with GLN189 and GLY143 showing an average bond length of 2.8 and 2.7 Å and an occupation percentage of 85.6 and 75.6%, respectively. Overall, these hydrogen bond analyses supported the finding of high stability for the erylosides B (226)-M pro complex compared to lopinavir-M pro .

To gain a more in-depth insight into the stability of ligand-M pro during the MD simulation, center-of-mass (CoM) distances were measured ( Figure 6 ). Interestingly, CoM distances were consistent for the erylosides B (226)-M pro complex compared to lopinavir-M pro , with average values of 5.6 vs. 5.9 Å, respectively. This suggests that erylosides B (226) bound more tightly to the M pro complex than lopinavir. 

To gain a more in-depth insight into the stability of ligand-M pro during t simulation, center-of-mass (CoM) distances were measured ( Figure 6 ). Interestingl distances were consistent for the erylosides B (226)-M pro complex compared to lop M pro , with average values of 5.6 vs. 5.9 Å, respectively. This suggests that erylo (226) bound more tightly to the M pro complex than lopinavir. 

To watch the structural stability of the erylosides B (226)-and lopina complexes, the root-mean-square deviation (RMSD) values of the backbone atom whole complex were estimated (Figure 7) . Unambiguously, the estimated RMSD for the investigated complexes remained below 0.20 nm over the 100 ns MD simu The erylosides B (226)-and lopinavir-M pro complexes reached the steady state in 10 ns MDs, and continued in an almost stationary state until the end of the simu The current results emphasized that the erylosides B (226) is tightly bonded and d impact the overall topology of SARS-CoV-2 M pro . for the investigated complexes remained below 0.20 nm over the 100 ns MD simu The erylosides B (226)-and lopinavir-M pro complexes reached the steady state in 10 ns MDs, and continued in an almost stationary state until the end of the simu The current results emphasized that the erylosides B (226) is tightly bonded and d impact the overall topology of SARS-CoV-2 M pro . 

Erylosides B (226) protein targets associated with severe acute respiratory syndrome diseases were predicted using a SwissTargetPrediction DisGeNET online tool. One hundred and seventeen genes were identified using Venn diagram comparison analysis. Commonly shared genes for erylosides B (226) were BCL2L1, IL2, PRKCA, and PRKCB ( Figure 8 ). Several studies reported that these four genes alter cytokine levels and immune functions in patients with COVID-19 infections and related conditions [26, 27] . Erylosides B (226) predicted gene targets were also analyzed via a STRING protein-protein interaction (PPI) network and visualized by Cytoscape 3.8.0. The top 20 scored genes for erylosides B (226) also included BCL2L1, IL2, PRKCA, and PRKCB (Table S3) . 

Erylosides B (226) protein targets associated with severe acute respiratory syndrome diseases were predicted using a SwissTargetPrediction DisGeNET online tool. One hundred and seventeen genes were identified using Venn diagram comparison analysis. Commonly shared genes for erylosides B (226) were BCL2L1, IL2, PRKCA, and PRKCB ( Figure 8 ). Several studies reported that these four genes alter cytokine levels and immune functions in patients with COVID-19 infections and related conditions [26, 27] . Erylosides B (226) predicted gene targets were also analyzed via a STRING protein-protein interaction (PPI) network and visualized by Cytoscape 3.8.0. The top 20 scored genes for erylosides B (226) also included BCL2L1, IL2, PRKCA, and PRKCB (Table S3 ). 

Toward a deep dissection and mining of the erylosides B (226) target-function interactions, PEA analysis and Boolean network modeling were conducted. A Voronoi treemap of erylosides B (226) top targeted pathway affected by the top 20 gene targets in response to erylosides B (226) in terms of SARS-CoV-2 infection was constructed ( Figure  9 ). Moreover, a Reactome hierarchy map was built to provide a genome-wide representation of the pathways influenced in response to erylosides B treatment ( Figure  S4 ). Remarkably, among the top 20 enriched pathways from the Reactome-PEA analyses and signal transduction, immune system and homeostasis signaling pathways were found to be the major pathways targeted by erylosides B (226), with a high significance (FDR < 0.00001%) (Table S4 ). (Figure 9 ). Moreover, a Reactome hierarchy map was built to provide a genome-wide representation of the pathways influenced in response to erylosides B treatment ( Figure S4 ). Remarkably, among the top 20 enriched pathways from the Reactome-PEA analyses and signal transduction, immune system and homeostasis signaling pathways were found to be the major pathways targeted by erylosides B (226), with a high significance (FDR < 0.00001%) (Table S4) . GPCRs are a massive family of surface transmembrane receptors in a human cell capable of responding to numerous stimulants and promoting a cascade of cellular functions. According to the GPCRs' pharmacological attributes, they are divided into four prominent families: the rhodopsin-like receptors family, the secretin receptors family, the metabotropic glutamate/pheromone receptors family, and the frizzled receptors family [28] . Remarkably, the Reactome-PEA analyses results revealed that almost all GPCR prominent families were significantly influenced by erylosides B (Figure 10 ). Many recent studies found that SARS-CoV-2 may compromise the GPCR signaling pathway and contribute to pulmonary edema's pathophysiology [29] . Additionally, recent reports speculated that SARS-CoV-2 might hijack the GPCR signaling pathway-mimicking a previous case report in cholera toxin, ADP-ribosylation Gαs-to activate CFTR and switch on Clsecretion, eventually leading to dysregulation of lung ion and fluids transmission, which ultimately led to lung edema in COVID-19 patients [30] . Interestingly, under signal transduction, it was found that the G-protein coupled receptor (GPCR) signaling pathway was the most enriched pathway influenced by erylosides B (226) treatment within the human genome. Mining of the Reactome-PEA analysis results emphasized that a set of 15 genes (OPRK1, HTR2C, PRKCA, HTR2A, ADRA1D, PRKCB, DRD2, DRD3, ADRA2A, EDNRB, ADRA1A, ADRA2B, ADRA1B, ADRA2C, and F2) were significantly modulated as biological targets to erylosides B (226) as potent SARS-CoV-2 inhibitor. Additionally, the Reactome-PEA results revealed that these genes were found to interact with other genes/interactors, including P16220, P52292, P00533, P12931, P07550, P35414, P01019, P30559, P14416, and O15354.

GPCRs are a massive family of surface transmembrane receptors in a human cell capable of responding to numerous stimulants and promoting a cascade of cellular functions. According to the GPCRs' pharmacological attributes, they are divided into four prominent families: the rhodopsin-like receptors family, the secretin receptors family, the metabotropic glutamate/pheromone receptors family, and the frizzled receptors family [28] . Remarkably, the Reactome-PEA analyses results revealed that almost all GPCR prominent families were significantly influenced by erylosides B (Figure 10 ). Many recent studies found that SARS-CoV-2 may compromise the GPCR signaling pathway and contribute to pulmonary edema's pathophysiology [29] . Additionally, recent reports speculated that SARS-CoV-2 might hijack the GPCR signaling pathway-mimicking a previous case report in cholera toxin, ADP-ribosylation Gαs-to activate CFTR and switch on Clsecretion, eventually leading to dysregulation of lung ion and fluids transmission, which ultimately led to lung edema in COVID-19 patients [30] . studies found that SARS-CoV-2 may compromise the GPCR signaling pathway and contribute to pulmonary edema's pathophysiology [29] . Additionally, recent reports speculated that SARS-CoV-2 might hijack the GPCR signaling pathway-mimicking a previous case report in cholera toxin, ADP-ribosylation Gαs-to activate CFTR and switch on Clsecretion, eventually leading to dysregulation of lung ion and fluids transmission, which ultimately led to lung edema in COVID-19 patients [30] . 

The resolved three-dimensional (3D) structure of SARS-CoV-2 main protease (M pro ) with a resolution of 2.16 Å (PDB code: 6LU7 [31] ) complexed with peptidomimetic inhibitor (N3) was downloaded and used as a template for all molecular docking and MD calculations. The protein structure was prepared by removing all heteroatoms, crystallographic waters, and ions, conserving only the amino acid residues. The protonation states of amino acid residues of M pro were assigned using H++ web server, and all missing hydrogen atoms were inserted [32] . The pKa values of M pro amino acid residues were estimated under physical conditions of salinity = 0.15, internal dielectric constant = 10, pH = 7, and external dielectric constant = 80.

The chemical structures of the 227 scrutinized marine natural products (MNPs) described in literature as Red-Sea terpenes were obtained in SDF format from the PubChem database (https://pubchem.ncbi.nlm.nih.gov). All MNPs were processed using Omega2 software to generate the three-dimensional (3D) structures of the studied molecules [33, 34] . Merck Molecular Force Field 94 (MMFF94S), implemented inside SZYBKI software [35, 36] , was applied to minimize the geometry of the investigated MNPs. The two-dimensional (2D) chemical structures of the scrutinized molecules are presented in Table S1 .

All molecular docking calculations were executed using AutoDock4.2.6 software [37] . The pdbqt file of SARS-CoV-2 main protease (M pro ) was prepared on the basis of the AutoDock protocol [38] . The maximum number of energy evaluations (eval) and the genetic-algorithm number (GA) were adjusted to 25,000,000 and 250, respectively. All other docking parameters were conserved at default. The grid box dimensions were adjusted to 60 Å × 60 Å × 60 Å to cover the SARS-CoV-2 M pro binding pocket. Moreover, the grid spacing value was set to 0.375 Å. The coordinates of the grid center were located at −13.069, 9.740, and 68.490 in the x, y, and z directions, respectively. The atomic partial charges of the investigated MNPs were calculated using the Gasteiger method [39] . The predicted binding poses for each investigated MNP were processed by the built-in clustering analysis (1.0 Å RMSD tolerance), and the conformation with the lowest energy within the most massive cluster was picked out as a representative pose.

AMBER16 software was employed to perform all MD simulations for the most potent MNPs complexed with SARS-CoV-2 M pro [40] . M pro and the investigated MNPs were determined using AMBER force field 14SB [41] and General AMBER force field (GAFF2) [42] , respectively. In the present study, implicit-solvent and explicit-solvent MD simulations were performed. In the implicit-solvent MD simulations, the atomic partial charges of the investigated MNPs were evaluated using the AM1-BCC method [43] . No periodic boundary conditions and no cut-off were applied for nonbonded interactions (specifically, a cutoff value of 999 Å was employed). Furthermore, the solvent influence was estimated using igb = 1 solvent model [44] . Energy minimization was first carried out on the docked MNP-M pro complexes for 500 steps, and the minimized complexes were then gently heated from 0 K to 300 K over 10 ps under NVT condition using a Langevin thermostat. Finally, the production stage was executed over 250 ps, and snapshots were assembled every 1 ps, giving 250 snapshots. In the present study, all implicit-solvent MD simulations were conducted using the CPU version of pmemd (pmemd.MPI) implemented inside AMBER16.

In explicit-solvent MD simulations, the atomic partial charges of the investigated MNPs were estimated using the restrained electrostatic potential (RESP) approach at the HF/6-31G* level with the help of Gaussian09 software [45, 46] . The TIP3P water model with periodic boundary conditions was utilized to solvate the MNP-M pro complexes in a cubic water box with a minimum distance to the box edge of 15 Å. Energy minimizations for 5000 steps on the solvated MNP-M pro complexes were performed using combined steepest and conjugate gradient methods. The minimized complexes were then smoothly heated from 0 to 300 K over 50 ps and a restraint of 10 kcal mol −1 Å −1 was applied on the main protease. The complexes were then equilibrated under NPT conditions for 1 ns. The production stage was then performed for each investigated M pro -inhibitor complex over simulation times of 10 ns, 50 ns, and 100 ns. The particle mesh Ewald (PME) method was applied to estimate the long-range electrostatic forces and energies. Cut-off for coulombic interactions with a value of 12 Å was considered [47] . To maintain the temperature at 298 K, Langevin dynamics with a gamma_ln collision frequency set to 1.0 was utilized. A Berendsen barostat with a pressure relaxation time of 2 ps was applied for the pressure control [48] . To constrain all bonds including hydrogen atoms, a SHAKE algorithm with a time step of 2 fs was utilized [49] . For binding energy calculations and post-dynamics analyses, energy values and coordinates were collected every 10 ps over the production stage. All explicit-solvent MD simulations were executed utilizing the GPU version of pmemd (pmemd.cuda) implemented inside AMBER16. All in silico calculations involving quantum mechanics, molecular docking calculations, and MD simulations were carried out on the CompChem GPU/CPU cluster (hpc.compchem.net). BIOVIA DS Visualize 2020 was utilized to generate all molecular graphics [50] .

The molecular mechanics/generalized Born surface area (MM/GBSA) approach [51] was applied to evaluate the binding free energies of the most potent MNPs complexed with SARS-CoV-2 M pro . To appoint the polar solvation energy, the modified GB model proposed by Onufriev [52] (igb = 2) was employed. The uncorrelated snapshots collected throughout the MD simulations were subjected to MM/GBSA (∆G binding ) energy calculations, which estimated as follows:

where the energy term (G) is evaluated as:

E ele and E vdw and are electrostatic and van der Waals energies, respectively. G SA is the nonpolar contribution to the solvation-free energy from the solvent-accessible surface area (SASA). G GB is the electrostatic solvation free energy calculated from the generalized Born equation. A single-trajectory approach was applied, in which the coordinates of each MNP-M pro , M pro , and MNP were acquired from a single trajectory. Because of the expensive computational request and low prediction thoroughness in most cases, entropy calculations were ignored [53, 54] .

According to the structural resemblance of recognized inhibitor-target integrations, target forecasting of the most auspicious marine natural products was carried out utilizing the online web-based tools of SwissTargetPredicition (http://www.swisstargetprediction. ch). Moreover, the DisGeNET online database (https://www.disgenet.org) was applied to gather the accessible database for Severe Acute Respiratory Syndrome (SARS) diseases. The InteractiVenn online tool was employed to design Venn diagram [55] . A functional database of STRING for top predicted targets was applied to generate protein-protein interaction (PPI) network [56] . Cytoscape 3.8.2 was applied to examine target-function relation according to the network topological [57] . Furthermore, to explore all probable target-function relations for the top 20 SARS disease-related genes based on their network profiles, pathway enrichment analysis was achieved using Cytoscape 3.8.2 [57] . Lastly, the ReactomeFIViz tool was integrated for modeling and visualizing all potential drug-target interactions [58] .

In the context of the COVID-19 pandemic, it is apparent that zoonotic diseases pose a huge risk to public health and economic stability. A lack of targeted treatments has spurred the exploration of novel drug leads for which computational approaches provide a relatively quick and cost-effective approach. Herein, marine natural products isolated from the Red Sea were screened in silico as potential M pro inhibitors. On the basis of molecular docking calculations, MD simulations, and molecular mechanics/generalized born surface area binding energy calculations, erylosides B (226) demonstrated a favorable binding affinity with ∆G binding < −51.0 kcal/mol with M pro . The stability of 226 complexed with M pro was confirmed using energetic and structural analyses over the 100 ns MD simulation. Ultimately, we hypothesize that the utilization of erylosides B may represent a potential therapeutic approach against COVID-19 due to its useful capability to modulate the GPCR signaling pathway in COVID-19 patients. Experimental validation is expected to further identify the suitability of natural products such as erylosides B to serve as a SARS-CoV-2 inhibitor.

Coronaviridae and SARS-associated coronavirus strain HSR1

Structural and biochemical characterization of the nsp12-nsp7-nsp8 core polymerase complex from SARS-CoV-2

Coronavirus disease 2019 (COVID-19): A literature review

Structural and functional properties of SARS-CoV-2 spike protein: Potential antivirus drug development for COVID-19

Viral targets for vaccines against COVID-19

SARS-CoV-2 evolution during treatment of chronic infection

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved alpha-ketoamide inhibitors

In silico prediction of potential inhibitors for the main protease of SARS-CoV-2 using molecular docking and dynamics simulation based drug-repurposing

In silico evaluation of prospective anti-COVID-19 drug candidates as potential SARS-CoV-2 main protease inhibitors

In silico drug discovery of major metabolites from spices as SARS-CoV-2 main protease inhibitors

In-silico drug repurposing and molecular dynamics puzzled out potential SARS-CoV-2 main protease inhibitors

Natural-like products as potential SARS-CoV-2 M(pro) inhibitors: In-silico drug discovery

Rutin and flavone analogs as prospective SARS-CoV-2 main protease inhibitors: In silico drug discovery study

Flavonoids: Promising natural compounds against viral infections

Potential bioactive glycosylated flavonoids as SARS-CoV-2 main protease inhibitors: A molecular docking and simulation studies

Flavonoids with inhibitory activity against SARS-CoV-2 3CLpro

Marine pharmacology in 2003-4: Marine compounds with anthelmintic antibacterial, anticoagulant, antifungal, anti-inflammatory, antimalarial, antiplatelet, antiprotozoal, antituberculosis, and antiviral activities; affecting the cardiovascular, immune and nervous systems, and other miscellaneous mechanisms of action

Marine pharmacology in 2000: Marine compounds with antibacterial, anticoagulant, antifungal, anti-inflammatory, antimalarial, antiplatelet, antituberculosis, and antiviral activities; affecting the cardiovascular, immune, and nervous systems and other miscellaneous mechanisms of action

Molecular architecture and biomedical leads of terpenes from red sea marine invertebrates

ENCORE: The effect of nutrient enrichment on coral reefs. Synthesis of results and conclusions

Isolation, antimicrobial and antitumor activities of a new polyhydroxysteroid and a new diterpenoid from the soft coral Xenia umbellata

Bioactive hydroperoxyl cembranoids from the Red Sea soft coral Sarcophyton glaucum

Lopinavir-ritonavir in patients admitted to hospital with COVID-19 (RECOVERY): A randomised, controlled, open-label, platform trial

Role of molecular dynamics and related methods in drug discovery

Molecular dynamics simulations in drug design

Altered cytokine levels and immune responses in patients with SARS-CoV-2 infection and related conditions

Network-based drug repurposing for novel coronavirus 2019-nCoV/SARS-CoV-2

Heterotrimeric G protein activation by G-protein-coupled receptors

Subcellular hot spots of GPCR signaling promote vascular inflammation

SARS-CoV-2 may hijack GPCR signaling pathways to dysregulate lung ion and fluid transport

Structure of M(pro) from SARS-CoV-2 and discovery of its inhibitors

H++: A server for estimating pKas and adding missing hydrogens to macromolecules

Conformer generation with OMEGA: Algorithm and validation using high quality structures from the Protein Databank and Cambridge Structural Database

OMEGA 2.5.1.4; OpenEye Scientific Software

MMFF94s option for energy minimization studies

AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility

Computational protein-ligand docking and virtual drug screening with the AutoDock suite

Iterative partial equalization of orbital electronegativity-A rapid access to atomic charges

Improving the accuracy of protein side chain and backbone parameters from ff99SB

Development and testing of a general amber force field

Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation

Implicit solvent models

A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges-The RESP model

Particle mesh Ewald: AnN·log(N) method for Ewald sums in large systems

Molecular-dynamics with coupling to an external bath

Settle-An Analytical Version of the Shake and Rattle Algorithm for Rigid Water Models

Combined molecular mechanical and continuum solvent approach (MM-PBSA/GBSA) to predict ligand binding

Exploring protein native states and large-scale conformational changes with a modified generalized born model

Assessing the performance of the molecular mechanics/Poisson Boltzmann surface area and molecular mechanics/generalized Born surface area methods. II. The accuracy of ranking poses generated from docking

End-Point Binding Free Energy Calculation with MM/PBSA and MM/GBSA: Strategies and Applications in Drug Design

InteractiVenn: A web-based tool for the analysis of sets through Venn diagrams

Anti-colorectal cancer targets of resveratrol and biological molecular mechanism: Analyses of network pharmacology, human and experimental data

Cytoscape: A software environment for integrated models of biomolecular interaction networks

Visualization of drug target interactions in the contexts of pathways and networks with ReactomeFIViz

The following are available online. Figure S1 : 2D representations of interactions of lopinavir and the top 27 potent marine natural products (MNPs) with the proximal amino acid residues of SARS-CoV-2 main protease (M pro ), Figure S2 : 3D representations of predicted binding modes of (i) 190, (ii) 178, (iii) 226 and (iv) lopinavir towards SARS-CoV-2 main protease (M pro ), Figure S3 : 3D representations of binding modes of (i) erylosides B (226)-and (ii) lopinavir-M pro complexes according to an average structure over a 100 ns MD simulation, Figure S4 : a genome-wide Reactome hierarchy map of the pathways influenced by the top 20 gene targets in response to erylosides B (226) in term of SARS-CoV-2 infection, Table S1 : evaluated docking score (in kcal/mol) for lopinavir and all investigated marine natural products (MNPs) against SARS-CoV-2 main protease (M pro ), Table S2 : computed Autodock and MM/GBSA binding energies (in kcal/mol) for the top 27 potent marine natural products (MNPs) against SARS-CoV-2 main protease (M pro ) over 250 ps implicit solvent MD simulations, Table S3 : network topological analysis for the predicted targets for erylosides B (226), Table S4: