1. Background
RNA viruses are a diverse group of human pathogens capable of triggering rapid seasonal outbreaks, posing serious public health challenges. COVID-19, caused by SARS-CoV-2, has recently evolved into a global pandemic with millions of deaths (1-3). Although primarily a respiratory illness, the disease also leads to various extrapulmonary complications due to the widespread presence of angiotensin-converting enzyme 2 (ACE2) receptors in different organs. These include thrombotic events, hyperglycemia, cardiac arrhythmias, acute coronary syndromes, kidney injury, gastrointestinal and liver dysfunctions, as well as neurological, ocular, and dermatological symptoms (4).
SARS-CoV-2 has a large positive-sense RNA genome comprising 14 open reading frames (ORFs). The 5′ ORFs (orf1a and orf1ab) encode polyproteins pp1a and pp1ab, which are cleaved into 16 nonstructural proteins (nsp1 - 16), while the 3′ region encodes structural proteins: Spike (S), envelope (E), membrane (M), and nucleocapsid (N) (5). Proteolytic cleavage of the polyproteins is catalyzed by two key enzymes: Three-chymotrypsin-like protease (3CLpro) and papain-like protease (PLpro) (6). Among them, 3CLpro, the main protease, consists of three domains — β-barrel domains I and II, and α-helical domain III — with the catalytic dyad HIS41 and CYS145 located between domains I and II (7).
Given the unique catalytic mechanism of 3CLpro, distinct from human proteases, its inhibition can effectively block viral replication, offering a promising target for specific antiviral therapies (8). Importantly, the high conservation of viral protease active sites reduces the likelihood of resistance development compared to vaccines (9). Although traditional drug discovery is often slow and expensive, natural compounds — especially those amenable to large-scale production — are gaining attention for their antiviral potential (10).
Marine organisms, which inhabit the majority of Earth’s surface, produce diverse secondary metabolites in response to environmental stressors. These compounds possess distinctive structural and functional characteristics, setting them apart from terrestrial natural products, and have shown promising antiviral activities. Given the rapid mutation and global spread of SARS-CoV-2, marine-derived agents represent a valuable resource for novel antiviral development (11).
Drug discovery is inherently complex, requiring interdisciplinary strategies to enhance efficiency and reduce cost. In this context, computational tools such as molecular docking and molecular dynamics (MD) simulations play a crucial role. Molecular docking predicts the interaction modes and affinities of small molecules with target proteins, while MD simulations provide insights into the dynamic behavior and conformational changes of protein-ligand complexes over time (12-14).
2. Objectives
In this study, we utilize these computational tools to explore the interactions between marine-derived compounds and the SARS-CoV-2 3CLpro protease. Through molecular docking and subsequent MD simulations, we identify promising inhibitors that may be developed into effective antiviral agents. The pharmacological properties of the most promising compounds are also evaluated to guide future experimental research.
3. Methods
3.1. Preparation of the Receptor and Ligands
The three-dimensional (3D) crystal structure of the SARS-CoV-2 main protease (3CLpro) was retrieved from the RCSB Protein Data Bank (PDB; ID: 7BAJ). The chemical structures of the selected marine-derived compounds (15, 16) were first sketched in two-dimensional (2D) format using the ACD/Labs software. These structures were then converted into 3D conformations and subjected to energy minimization using the steepest descent algorithm in the Avogadro software.
3.2. Molecular Docking Studies
Molecular docking was employed to predict the binding affinity and orientation of ligands within the active site of 3CLpro (17, 18). Docking simulations were performed using AutoDock Vina. Prior to docking, all ligand torsions were set as rotatable, and receptor preparation, including the addition of polar hydrogens, assignment of Gasteiger charges, and atom typing according to the AutoDock force field, was carried out using MGLTools. The docking grid was centered on the protease's catalytic dyad (His41 and Cys145), defining the active site. Each of the 80 marine compounds was subjected to 20 docking runs. Ligands exhibiting binding energies lower than -7 kcal/mol (Table 1) and interacting with the active site residues were selected for further evaluation of pharmacokinetic properties and drug-likeness using the SwissADME web server (http://www.swissadme.ch). Based on these assessments, five promising candidates with proper interaction energies and pharmacokinetic properties were chosen for MD simulations.
| Compound | Lowest Binding Energy (kcal/mol) |
|---|---|
| 3-methoxydebromoaplysiatoxin | -8.2 |
| Alterporriol Q | -7.9 |
| Arisugacin A | -7.7 |
| Asebotin | -7.5 |
| Aspernolide A | -7.4 |
| Asperterrestide A | -7.4 |
| Baculiferins k | -8.4 |
| BaculiferinsC | -8.4 |
| BaculiferinsE | -7.5 |
| BaculiferinsF | -8 |
| BaculiferinsH | -7.2 |
| BaculiferinsL | -8.1 |
| BaculiferinsM | -8.3 |
| BaculiferinsN | -8.5 |
| briacavatolides E | -7.6 |
| Cadiolide B | -8.3 |
| Cordyol C | -7 |
| Diekol | -8 |
| DurumolideJ | -7.8 |
| Ehrenbergol C | -7 |
| Gyrosanol A | -7.9 |
| Halistanolsulfate | -7.3 |
| HipostroneN | -7.4 |
| Isobutyrolactone II | -7.4 |
| Khayanolides | -7.2 |
| Manoalide | -7.2 |
| Manzamine A | -7.5 |
| MirabamidesH | -7.4 |
| Mollamide E | -7.1 |
| Mollamide F | -7.5 |
| Molleurea A | -7.1 |
| Myticin | -7 |
| Norlichexanthone-3,6,8-Trihydroxy-1-methylxanthone | -7.4 |
| Nortopsentins | -7.6 |
| Penipanoid C | -7.4 |
| Phlorofucofuroeckol | -7.6 |
| Polyhydroxylated steroids | -7.5 |
| Polyhydroxylated sterol | -7.7 |
| Prunolide A | -7.5 |
| Quercetin-3-O-β-D-xylopyranoside | -8.6 |
| Rubrolide S | -7.5 |
| Sorbicatechols A | -7.2 |
| Stachybotrins D | -7.9 |
| ThalassiolinD | -8 |
| Thalassodenrone | -7.4 |
| Triterpenoids | -7.5 |
Binding Energies of 45 Selected Marine-derived Compounds Against SARS-CoV-2 Main Protease (3CLpro) from Molecular Docking Simulations
3.3. Molecular Dynamic Simulation
The MD simulations were performed using the GROMACS 2018 package (19) to assess the structural stability and dynamic behavior of five ligand-3CLpro complexes. Ligand topologies were prepared via the Antechamber software, and the protein topology was generated using the AMBER99SB-ILDN force field (20). Each complex was solvated in a cubic box with TIP3P water molecules, and counterions were added for charge neutrality. Energy minimization was achieved through the steepest descent method until atomic forces dropped below 10 kJ/mol/nm (21). Periodic boundary conditions were applied in all spatial dimensions, and equilibration was performed under both NVT and NPT ensembles at 310 K and 1 bar using the V-rescale thermostat (22) and Martonak et al. (23), respectively. Short-range electrostatic and van der Waals interactions were calculated with a 1.0 nm cutoff, and bond constraints were applied using the LINCS algorithm (24). Simulations were run for 100 ns using the leap-frog integrator (25). Post-simulation analyses, including structural evaluations and trajectory visualizations, were carried out using VMD, and 2D interaction diagrams were generated with LigPlot+ (26).
4. Results and Discussion
4.1. Molecular Docking Studies
Molecular docking was utilized to assess the binding affinity and orientation of 80 marine-derived natural compounds against the SARS-CoV-2 main protease (3CLpro). Compounds with binding energies below -7 kcal/mol were considered promising, resulting in the selection of 45 candidates for further pharmacokinetic and drug-likeness evaluation via the SwissADME server. Docking scores are listed in Table 1. Based on the docking and ADME results, five top-performing ligands (Appendices 1 and 2 in Supplementary File) — 3-methoxydebromoaplysiatoxin, Aspernolide A, Ehrenbergol C, Isobutyrolactone II, and Penipanoid C — were selected for detailed interaction analysis and subsequent molecular dynamics simulations.
Docking interaction visualization (Figure 1) revealed the detailed binding of 3-methoxydebromoaplysiatoxin to active site residues, including the catalytic dyad His41 and Cys145, through multiple hydrogen bonds and hydrophobic contacts. Aspernolide A formed critical interactions with domain I residues and Glu166, engaging the catalytic site and suggesting enzyme interference. Ehrenbergol C showed stable hydrogen bonding with residues in domains I and II, including Cys145, implying partial inhibition of substrate access. Isobutyrolactone II interacted strongly with both catalytic residues and several nearby residues, indicating potential enzymatic blockage. Penipanoid C established hydrogen and van der Waals interactions with residues essential to catalysis, possibly disrupting substrate cleavage. Overall, the docking results indicate that these marine compounds have a strong affinity for the 3CLpro active site, justifying further investigation through MD simulations.
Two-dimensional (2D) and three-dimensional (3D) visualizations of selected marine compounds docked with the active site of SARS-CoV-2 3-chymotrypsin-like protease (3CLpro): A, 3-Methoxydebromoaplysiatoxin; B, Aspernolide A; C, Ehrenbergol C; D, Isobutyrolactone II; and E, Penipanoid C.
4.2. Molecular Dynamics Simulations
Molecular dynamics simulations were performed using 100 ns all-atom trajectories to study the conformational stability and interaction patterns of selected marine-derived compounds with the SARS-CoV-2 main protease (3CLpro). Analyses included assessments of changes in binding mode and persistence (Appendix 3 in Supplementary File), interaction dynamics, and structural changes over time. For 3-Methoxydebromoaplysiatoxin, hydrogen bonds were reduced to two with Thr26 and Met49 in domain I, and the compound retained interaction only with Cys145, indicating a slight shift from its original binding pose. Aspernolide A preserved interactions with domain II residues (Met165, Val142, Arg188, Gln189) and maintained two hydrogen bonds with His41 and Glu166, although the total number of interactions declined. Ehrenbergol C lost several original contacts and retained only hydrophobic interactions with Glu166, Gln189, and Thr190 in domain II, suggesting weakened binding. Isobutyrolactone II established a strong hydrogen bond with Glu47 and retained interactions with domain I residues and the catalytic dyad His41 and Cys145, supporting its stable binding. Penipanoid C also showed reduced residue interactions post-simulation but preserved contacts with catalytic residues His41 and Cys145, formed a stable hydrogen bond with Glu166, and retained van der Waals interactions with Met49, Met165, and Gln189, indicating potential to inhibit enzymatic function. Collectively, these findings demonstrate that while marine-derived ligands undergo subtle conformational adjustments during MD simulations, they retain stable binding interactions with the active site of the protease. The persistent structural integrity of the complexes, coupled with sustained ligand-enzyme contact at the catalytic site, suggests a plausible mechanistic basis for their protease inhibitory activity.
Root mean square deviation (RMSD) analysis over 100 ns for protein backbone atoms revealed that the unbound protease initially increased in RMSD during the first 20 ns and then stabilized (Figure 2). The 3-Methoxydebromoaplysiatoxin complex showed its highest RMSD peak near 40 ns and stabilized around 0.36 nm thereafter. The Aspernolide A complex fluctuated in the first 20 ns, then stabilized around 0.33 nm, about equal to the free protein from 40 ns onward. The Isobutyrolactone II and Penipanoid C complexes both showed rising RMSD trends. Penipanoid C and Isobutyrolactone II had values exceeding 0.42 nm, especially between 70 - 100 ns, whereas Ehrenbergol C had the lowest average RMSD. These results suggest that different ligands exert distinct effects on the protease's conformational dynamics.
Root mean square deviation (RMSD) profiles over time for five protein-ligand complexes: A, 3-Methoxydebromoaplysiatoxin; B, Aspernolide A; C, Ehrenbergol C; D, Isobutyrolactone II; and E, Penipanoid C. The RMSD profile of the free protein is shown in pink for comparison.
Figure 3 illustrates how the binding of marine-derived ligands to the 3CLpro influences residue flexibility across its domains. In the 3-Methoxydebromoaplysiatoxin-3CLpro complex, root mean square fluctuation (RMSF) decrements were observed at residues in domains I (41 - 45, 46 - 52, 59 - 70, 84), II (139 - 142, 180 - 187), and III (188 - 193). Along with increased fluctuations in residues 7 - 15, 91 - 100 (domain I), 101 - 111, 118 - 124, 128 - 137 (domain II), and 223 - 234, 236 - 258, 275 - 306 (domain III), these suggested variable changes in the flexibility of domains I and II, and increased flexibility in domain III.
Root mean square fluctuation (RMSF) profiles of Cα atoms for five protein–ligand complexes: A, 3-Methoxydebromoaplysiatoxin; B, Aspernolide A; C, Ehrenbergol C; D, Isobutyrolactone II; and E, Penipanoid C. The RMSF profile of the free protein is shown in pink for comparison.
For Aspernolide A, RMSF increased at residues in domains I (19 - 24, 55 - 65, 92 - 98), II (105 - 115, 116 - 134, 137 - 141, 158 - 165), and III (188 - 191, 274 - 279), while decreases occurred at residues 7 - 17, 50 - 52, 72, 76 - 87 (domain I), 177 - 183, 185 - 190 (domain II), and across most of domain III. The Ehrenbergol C complex showed RMSF decreases at residues 43 - 50, 63 - 65, 72, 83 - 86 (domain I), 140 - 142, 169 - 174, 177 - 181 (domain II), and 182 - 194, 224 - 226, 229 - 239, 252 - 263 (domain III). Increases were found at residues 106 - 116, 120 - 134, 155 (domain II) and 202 - 220, 245 - 252, 275 - 306 (domain III). Isobutyrolactone II caused increased RMSF at residues in domains I (67 - 70, 75 - 83), II (88 - 100, 102 - 116, 118 - 145, 153 - 171), and III (200 - 207, 213 - 227, 237 - 306), while decreases were noted at 46 - 50, 56, 61 - 63 (domain I), 183 - 196 (domain II), and 208 - 211 (domain III). Penipanoid C induced RMSF decreases at residues 14, 41 - 52, 54 - 59, 64 (domain I), 115, 119 - 135, 139 - 142, 168 - 182 (domain II), and 183 - 196, 209, 231 - 235 (domain III), with increases at 71, 75 - 80, 91 - 98 (domain I), 102 - 112 (domain II), and 203 - 306 (domain III). In conclusion, ligand binding significantly affects the flexibility of key regions in 3CLpro involved in catalysis and substrate binding, supporting the inhibitory potential of these marine-derived compounds.
The radius of gyration (Rg) analysis measures the overall changes in the diameter of a protein. Changes in Rg indicate conformational variations that could affect substrate binding to the catalytic pocket. The Rg profiles of the free protein and its complexes with various marine compounds are shown in Appendix 4. The free protein showed some compactness with fluctuations around an average of 2.16 nm, showing a downward trend during the simulation. In contrast, complexes with all ligands prevent the final Rg reductions in ligand-bonded protein when compared with the free protein. The highest fluctuations were seen in complexes with Ehrenbergol C, Isobutyrolactone II, Aspernolide A, and 3-Methoxydebromoaplysiatoxin, suggesting more structural instability and potential enzyme dysfunction.
Solvent-accessible surface area (SASA) analysis, shown in Appendix 5, revealed a small increase in the SASA value in the Isobutyrolactone II complex, indicating a structural change of the protein in the presence of this ligand. As shown in the figure, in other complexes, the final SASA value is in the range of the free protein with some differences in diagram patterns. The most severe fluctuations in SASA value during the simulation are observed in the protein complex with Penipanoid C and 3-Methoxydebromoaplysiatoxin.
Principal component analysis (PCA) extracted the protein’s movement patterns, as shown in Figure 4. The free protein exhibited a value range of -7 to 8 nm in PC1. All marine compounds altered the range and pattern of conformational values in PCs, indicating their potential to destabilize the structure. Notably, Penipanoid C caused the largest increase in protein mobility, consistent with the Rg results, showing significant fluctuations. In contrast, complexes with Isobutyrolactone II, Aspernolide A, and Ehrenbergol C showed reduced protein mobility, with the most severe restriction observed in Isobutyrolactone II. These findings suggest that these marine compounds could cause conformational changes and movement patterns of the protease, hindering the protein's activity.
Principal component analysis (PCA) patterns of: A, free protein, and protein complexes with B, 3-Methoxydebromoaplysiatoxin; C, Aspernolide A; D, Ehrenbergol C; E, Isobutyrolactone II; and F, Penipanoid C.
The Define Secondary Structure of Proteins (DSSP) analysis (Appendix 6 in Supplementary File) demonstrated changes in the protein's secondary structures resulting from ligand binding. For example, 3-methoxydebromoaplysiatoxin destabilized an alpha-helix between residues 50 - 60 and removed a turn between 170 - 180. Aspernolide A disrupted beta - sheets between residues 90 - 130, while Ehrenbergol C removed a turn at residues 40 - 50. Isobutyrolactone II induced an alpha-helix formation between residues 10 - 20 and removed a turn at 40 - 50. Penipanoid C promoted alpha-helix formation at residues 10 - 20 and 50 - 60, and coils at 190 - 200. These alterations suggest that ligand binding impacts the protein’s secondary structure and may affect its enzymatic activity.
Pharmacokinetic analysis of the marine compounds (Appendices 1 and 2 in Supplementary File), performed using the Swiss ADME server, revealed that all compounds exhibited good drug-likeness. Notably, Ehrenbergol C and Isobutyrolactone II did not interact with hepatic metabolism enzymes, suggesting a lower likelihood of drug interactions. Aspernolide A and Penipanoid C interacted with CYP2C9, CYP3A4, and CYP1A2, potentially interfering with anti-inflammatory drug metabolism. All five compounds exhibited high digestive absorption and the ability to cross the blood-brain barrier, suggesting efficient distribution and potential effectiveness in reaching viral enzymes.
4.3. Conclusions
This study computationally assessed 80 marine-derived compounds for their inhibitory potential against the SARS-CoV-2 main protease (3CLpro). Based on binding affinity and ADMET analysis, five lead compounds — 3-methoxydebromoaplysiatoxin, Aspernolide A, Ehrenbergol C, Isobutyrolactone II, and Penipanoid C — were selected for MD simulations.
The MD results showed that complexes with Isobutyrolactone II, Ehrenbergol C, and 3-methoxydebromoaplysiatoxin exhibited the highest RMSD values, with Isobutyrolactone II causing the greatest fluctuations, confirmed by both RMSD and residue-level analysis. The Rg analysis indicated notable diagram fluctuations in 3-methoxydebromoaplysiatoxin, Ehrenbergol C, and Isobutyrolactone II complexes, while others remained structurally more stable. The results of PCA revealed altered protein motion in all complexes, with the most significant limitation in protein movement observed for Isobutyrolactone II, Ehrenbergol C, and Aspernolide A, likely affecting enzymatic function. SwissADME pharmacokinetic analysis identified Isobutyrolactone II and Penipanoid C as having the most favorable pharmacokinetic profiles due to favorable absorption and metabolic characteristics. Integrative analysis of computational data nominates Isobutyrolactone II and Aspernolide A as high-priority candidates for SARS-CoV-2 antiviral development. Their persistent interaction with the viral protease active site, alongside their negative effects on the normal dynamics of the protein, underscores their inhibitory potential, necessitating further investigation in biological systems to confirm therapeutic utility.
