Identification of potential Indonesian marine invertebrate bioactive compounds as TMPRSS2 and SARS-CoV-2 Omicron spike protein inhibitors through computational screening

2.1 Collection of macromolecules

The human TMPRSS2 crystal structure complex with nafamostat as an inhibitor (PDB ID: 7MEQ) and SARS-CoV-2 Omicron spike protein B.1.1.529 (PDB ID: 7QO7) three-dimensional structures were obtained from the Research Collaboratory for Structural Bioinformatics Protein Data Bank website (https://www.rcsb.org/). TMPRSS2 is a single-chain proenzyme (Fraser et al. 2022), while the SARS-CoV-2 Omicron spike protein structure is a trimer spike with monomers A, B, and C (Ni et al. 2021). Table 1 outlines the structural characteristics of TMPRSS2 and SARS-CoV-2 Omicron spike protein.

2.2 Preparation and optimization of macromolecules

Unwanted ligands and TMPRSS2 and SARS-CoV-2 Omicron spike chain C structures were removed from water molecules. The macromolecules were visually examined using the Chimera tool to ensure that the empty amino acid residues were filled. Modeller 10.3 was applied to create all of the missing gaps and loops (Fiser et al., 2000). The ModLoop web server (https://modbase.compbio.ucsf.edu/modloop) was used to model and enhance the loops. Autodock Tools 1.5.6 was utilized to incorporate polar hydrogens and Kollman united atom charges in the 3D structures of the proteins (Trott and Olson 2009).

2.3 Collection and preparation of ligands

Overall, 137 secondary metabolites of Indonesian marine invertebrates were utilized as test ligands (Izzati et al. 2021; Nurrachma et al. 2021). Camostat was used as a reference inhibitor of TMPRSS2, and mefloquine was utilized as a reference inhibitor of the spike protein.

MarvinSketch 20.14 was employed to sketch the test ligand molecules, and the three-dimensional structures were created by adding polar hydrogen. The 3D structure of the ligand molecule was created by adding hydrogen atoms using BIOVIA Discovery Studio Visualizer 2021 software. As reference ligands, the 3D structures of camostat and mefloquine were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Meanwhile, the nafamostat’s 3D structure was obtained by separating the nafamostat co-crystal bound to the TMPRSS2 protein.

2.4 Determination of the binding site of SARS-CoV-2 Omicron spike protein with TMPRSS2: Protein-Protein Interaction (PPI)

The structures of the SARS-CoV-2 Omicron spike protein and TMPRSS2 were docked by the HADDOCK 2.4 web server (https://wenmr.science.uu.nl/haddock2.4/) based on a previous study by Salleh and Deris (2022). The model with the best Z-HADDOCK score and the largest number of clusters was chosen to investigate the binding energy (Delta G) and dissociation constant (K_D) utilizing the PRODIGY web server (https://wenmr.science.uu.nl/ prodigy/). Predictions of surface contact hydrogen bonds between proteins were carried out via the PDBePISA web server (https://www.ebi.ac.uk/pdbe/pisa/).

2.5 Virtual screening

Virtual screening of the test ligands to TMPRSS2 and SARS-CoV-2 Omicron spike protein was performed using PyRx-Autodock Vina 0.9.8 (Trott and Olson 2009). The coordinates of the TMPRSS2 active site used were ×  = -9.755 Å, y = -6.247 Å, and z = 20.307 Å, with a grid box size of 68 × 68 × 68 pts and a grid spacing of 0.375 Å. According to the protein–protein interaction (PPI) data, the coordinates of the SARS-CoV-2 Omicron spike protein were based on the amino acid residues of the spike protein that were in contact to create hydrogen bonds with the TMPRSS2 protein. The virtual screening binding affinity energy scores were sorted, and the top 20 test ligand complexes were examined using the molecular docking method.

2.6 Molecular docking

The top 20 ligand–protein complexes identified through virtual screening were confirmed using molecular docking with the AutoDock Tools program. Validation was acceptable if the RMSD was less than 2 Å (Allen and Rizzo 2014). The number of Genetic Algorithm runs was set to 50, the population size was 250, the maximum number of evals was 2,500,000 (medium parameter), and the maximum number of generations was 27,000 (Gurung et al. 2022). The 10 docking-validated test ligands with the best binding affinity energy values were examined for binding affinity energy (Delta G), inhibition constants (Ki), and ligand interactions with amino acid residues. The BIOVIA Discovery Studio Visualizer tool was used to visualize the complexes’ test ligand–protein and reference ligand–protein interactions of amino acids.

2.7 ADMET prediction analysis

The ADMETlab 2.0 web server (https://admetmesh.scbdd.com) was used to analyze drug similarities and estimate absorption, distribution, metabolism, excretion, and toxicity profiles (Sharma et al. 2022). Lipinski’s rules, which consider molecular weight, the log P partition coefficient value, the number of hydrogen bond donors, and the number of hydrogen bond acceptors, were used to see if the compounds were similar to drugs. In addition to Lipinski’s rules, the absorption, distribution, metabolism, excretion, and toxicity (ADMET) profiles of bioactive compounds must fulfil the following criteria: topological polar surface area (TPSA) ≤ 75 Å², human intestinal absorption (HIA) ≤ 0.7, half-life (T_1/2) ≤ 0.8, carcinogenicity ≤ 0.3, and AMES toxicity ≤ 0.3 (Lazniewski et al. 2022). Molecular dynamics simulation was examined on the best test ligand that met all of the ADMET analysis requirements.

2.8 Molecular dynamic simulation

After selecting a ligand molecule based on ADMET profile analysis, molecular dynamic simulation using AMBER MD software was performed in multiple stages. First, macromolecules and ligands were prepared, and then the topologies and coordinates of macromolecules, ligands, and macromolecule-ligand complexes were established in a vacuum and in aqueous solvents. The next phase was energy minimization and equilibration. During the equilibration step, the volume was held constant, and the temperature was increased to 310 K. Three factors needed to be verified before production began. Temperature, density, and potential energy must all be held constant, with a temperature of 310 K and a density of 1 g/mL. Simulations for 100 ns were performed for the production to be generated with a stable balanced system (Sharma et al. 2022). The Cpptraj subprograms were used to assess the findings of the molecular dynamics simulation. Root mean square deviation (RMSD) and root mean square fluctuation (RMSF) were the parameters studied.

3.1 Protein-Protein interaction: determination binding site of SARS-CoV-2 Omicron spike protein

The SARS-CoV-2 Omicron spike protein was docked with TMPRSS2 to find where the spike protein binds. This information was then used in virtual screening with test ligands. When the SARS-CoV-2 spike protein binds to TMPRSS2, cleavage occurs in the S2 subunit (Ma et al. 2021), which activates the membrane fusion of the SARS-CoV-2 spike protein into the host cell (Evans and Liu 2021; Simmons et al. 2013; Zhu et al. 2021). The test ligand can prevent protease cleavage and thereby prevent SARS-CoV-2 from entering cells via endocytosis by preventing the interaction of the SARS-CoV-2 spike protein with TMPRSS2.

Protein-protein docking yielded 10 clusters with a total of 161 structures, accounting for 80% of the water-refined HADDOCK model. Cluster 1 had the best HADDOCK score of −69.3 ± 4.0, the largest cluster size of 38, and a Z-score of −0.8 (Supplementary File).acknow

The Z-score indicates the number of standard deviations from the cluster mean; the lower the value, the better. The cluster 1 model was then used to calculate binding energy (Delta G) and dissociation constant (KD) utilizing the PRODIGY web server (https://wenmr.science.uu.nl/prodigy/). The investigation yielded binding energies (Delta G) of −9.5 kcal/mol and dissociation constants (KD) of 1.1 × 10–7 M. The PDBePISA web server was used to predict surface contact hydrogen bonding between proteins. The surface contact hydrogen bonds formed between the SARS-CoV-2 Omicron spike protein and TMPRSS2 can be seen in Table 2.

The BIOVIA Discovery Studio Visualizer program was used to depict the cluster 1 model (Fig. 3). TMPRSS2 bonded only to a single binding site in the active residues 685 and 686 of the 10 clusters of the SARS-CoV-2 Omicron spike protein. The bonds between TMPRSS2 and spike protein were enhanced by the formation of many h-bonds, one of which was the bond produced by GLY 462 from TMPRSS2 and SER 686 from the spike proteins’ cleavage site 1. Although the cleavage site of the spike protein did not directly contact all TMPRSS2 binding sites, the amino acid residues of the spike protein engaged with the residues of amino acids of TMPRSS2 clamped to the TMPRSS2 active site. Once attached, this flexible loop can alter conformation, bringing cleavage sites 685 and 686 of the spike protein nearer to the catalytic active site of TMPRSS2. In the two-step mechanism, this results in the enzymatic cleavage of S1/S2 and S’ and the fusion of the host cell membrane (Saxena et al. 2020).

Close

Fig. 3

Result of TMPRSS2 protein docking (PDB ID: 7MEQ) with Omicron variant SARS-CoV-2 spike protein (PDB ID: 7QO7) at cleavage site 1 active residue 686 and 685 (red residues). The SARS-CoV-2 Omicron variant spike protein is represented by magenta, while TMPRSS2 is represented by blue.

Export to PPT

3.2 Docking validation of TMPRSS2 protein and SARS-CoV-2 Omicron spike protein

Docking method validation is acceptable if the RMSD value less than 2 Å (Allen and Rizzo 2014). The deviation distance from the co-crystal or native ligand bond location with the protein after redocking to the actual co-crystal bond position is defined as the RMSD. The TMPRSS2 protein (7MEQ) and the SARS-CoV-2 Omicron spike protein (7QO7) were validated by redocking the reference ligands.

The TMPRSS2 protein worked by redocking nafamostat, a native ligand for TMPRSS2. The parameter of the grid box was 68 × 68 × 68 pts, the grid spacing was 0.375 Å, and the coordinates corresponding to the bond site were ×  = -9.755 Å, y = -6.247 Å, and z = 20.307 Å. In the same way, mefloquine, a reference ligand, was docked to validate the SARS-CoV-2 Omicron spike protein. The parameter of the grid box was 76 × 76 × 76 pts, grid spacing was 0.375 Å, and the coordinates were ×  = 233.664 Å, y = 171.257 Å, and z = 197.310 Å. This grid box was set up based on the amino acid residues on the spike protein that interact with TMPRSS2 to form hydrogen bonds. The PyMOL tool was used to align nafamostat and mefloquine ligands before and after redocking, and the findings indicate that the RMSD values were 0.424 Å and 0.057 Å, respectively (Fig. 4A, B). Based on these findings, the docking method’s validation produces good results, is acceptable, and could be utilized for virtual ligand screening tests.

Close

Fig. 4

Alignment of nafamostat (A) and mefloquine (B) for validation of the TMPRSS2 protein and the SARS-CoV-2 Omicron spike protein, respectively. The light blue color indicates the position of the ligand before redocking, while the light green color shows its position after redocking.

Export to PPT

3.3 Molecular docking-based virtual screening and ADMET analysis

According to the virtual screening and docking validation results, each of the 10 compounds exhibited the highest binding affinity energy for the TMPRSS2 protein and the SARS-CoV-2 Omicron spike protein, as shown in Tables 3 and 4. Compounds with affinity energies ranging from −10.43 to −8.63 kcal/mol were discovered during the screening of test ligands with TMPRSS2 (Table 3). The bond affinity energy of the test ligand exceeds that of the nafamostat co-crystal, which is −6.52 kcal/mol. Camostat has a lower comparable ligand bonding affinity energy (-8.25 kcal/mol) than nafamostat. In contrast, the findings of screening test ligands with SARS-CoV-2 Omicron spike protein (Table 4) revealed that compounds had an affinity energy range of −9.68 to −8.70 kcal/mol. The binding affinity energy of the tested ligands was higher than the binding affinity energy of the reference ligand mefloquine, which was −6.34 kcal/mol. Molecular docking-based virtual screening revealed that the test ligands derived from the secondary metabolites of Indonesian marine invertebrates were more successful than the reference ligands at inhibiting TMPRSS2 and SARS-CoV-2 Omicron spike protein. However, more research must be done to ensure that the hit compounds work against COVID-19 in vitro and in vivo.

A novel drug’s pharmacokinetics and toxicity must be known before it can be developed. The results show that acanthomanzamine C, cortistatin G, cortistatin J, and cortistatin L met the ADMET analysis requirements (Table 5) and indicated feasible oral administration. TPSA (Topological Polar Surface Area) represents a molecule's total surface of polar atoms, while HIA (Human Intestinal Absorption) shows how well the drug is absorbed in the gastrointestinal tract. TPSA and HIA parameters are useful indicators to describe the bioaccessibility and bioavailability of a compound. Parameter T1/2 (half-life) is used to set the interval of drug administration. Carcinogenicity and AMES toxicity are important to know whether a compound has carcinogenic or mutagenic effects. The BIOVIA Discovery Studio Visualizer tool was used to depict the top two test ligands and their reference ligands with the greatest binding affinity energy to each TMPRSS2 (Fig. 5) and SARS-CoV-2 Omicron spike protein (Fig. 6). The details of the H-bond interactions and those with hydrophobic interactions can be seen in Tables 3 and 4.

Table 5 ADMET analysis results.

Close

Fig. 5

Interactions of the reference ligands (camostat and nafamostat) and the test ligands (acanthomanzamine C and cortistatin G) with TMPRSS2 amino acid residues.

Export to PPT

Close

Fig. 6

Interactions of the reference ligand (mefloquine) and the test ligands (acanthomanzamine C and cortistatin J) with SARS-CoV-2 Omicron spike protein amino acid residues.

Export to PPT

3.4 Molecular dynamics simulation studies

Acanthomanzamine C, the primary molecule in the complex with TMPRSS2 and SARS-CoV-2 Omicron spike protein, was subjected to molecular dynamics simulations to determine the extent of this compound’s interactions with TMPRSS2 and SARS-CoV-2 Omicron spike protein. The most important ways to evaluate protein dynamics are RMSD and RMSF analyses.

The dynamics simulation of the ligands nafamostat (co-crystal), camostat (reference ligand), and acanthomanzamine C in association with TMPRSS2 protein was examined (Fig. 7A). The combination of nafamostat and TMPRSS2 was shown to be stable within the 100-ns period with RMSD less than 3 Å. In comparison to the nafamostat, the camostat and acanthomanzamine C exhibited an increase in RMSD after 40 ns, after which the graph stabilized up to 100 ns at RMSD less than 4 Å. Fig. 7B depicts the RMSD mefloquine (reference ligand) and acanthomanzamine C complexed with the SARS-CoV-2 Omicron spike protein. Mefloquine and acanthomanzamine C showed slight fluctuations but stabilized within the first 50 ns, with small variations found in the remaining duration at RMSD less than 4 Å. The RMSF figure (Fig. 8A, B) shows a high degree of variation in the amino acid residues implicated in the loop region of TMPRSS2 and SARS-CoV-2 Omicron spike protein, with maximum degrees of fluctuation reaching 6 Å and 4 Å, respectively. A careful examination reveals that the residues implicated in the ligand’s binding have low RMSF values.

Close

Fig. 7

RMSD results for 100 ns: A) ligands in a complex with the TMPRSS2 and B) ligands in a complex with the SARS-CoV-2 Omicron spike protein.

Export to PPT

Close

Fig. 8

RMSF results for 100 ns: A) amino acid residues of the TMPRSS2 and B) amino acid residues of the SARS-CoV-2 Omicron spike protein.

Export to PPT