Marine Brominated Tyrosine Alkaloids as Promising Inhibitors of SARS-CoV-2

There have been more than 150 million confirmed cases of SARS-CoV-2 since the beginning of the pandemic in 2019. By June 2021, the mortality from such infections approached 3.9 million people. Despite the availability of a number of vaccines which provide protection against this virus, the evolution of new viral variants, inconsistent availability of the vaccine around the world, and vaccine hesitancy, in some countries, makes it unreasonable to rely on mass vaccination alone to combat this pandemic. Consequently, much effort is directed to identifying potential antiviral treatments. Marine brominated tyrosine alkaloids are recognized to have antiviral potential. We test here the antiviral capacity of fourteen marine brominated tyrosine alkaloids against five different target proteins from SARS-CoV-2, including main protease (Mpro) (PDB ID: 6lu7), spike glycoprotein (PDB ID: 6VYB), nucleocapsid phosphoprotein (PDB ID: 6VYO), membrane glycoprotein (PDB ID: 6M17), and non-structural protein 10 (nsp10) (PDB ID: 6W4H). These marine alkaloids, particularly the hexabrominated compound, fistularin-3, shows promising docking interactions with predicted binding affinities (S-score = −7.78, −7.65, −6.39, −6.28, −8.84 Kcal/mol) for the main protease (Mpro) (PDB ID: 6lu7), spike glycoprotein (PDB ID: 6VYB), nucleocapsid phosphoprotein (PDB ID: 6VYO), membrane glycoprotein (PDB ID: 6M17), and non-structural protein 10 (nsp10) (PDB ID: 6W4H), respectively, where it forms better interactions with the protein pockets than the native interaction. It also shows promising molecular dynamics, pharmacokinetics, and toxicity profiles. As such, further exploration of the antiviral properties of fistularin-3 against SARS-CoV-2 is merited.


Introduction
The 2019 novel coronavirus disease , caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has significantly impacted global health and economics [1]. The signs and symptoms of COVID-19 are grouped into three categories according to the severity of the infection and mortality: mild, severe, and critical. The majority of COVID-19 patients, 80%, experience mild symptoms and recover. Severe symptoms appear in 13.8% of cases and 6.1% become critically ill [2,3]. SARS-CoV-2 is the seventh coronavirus known to infect humans, but the only one, to date, that has caused a pandemic [4,5]. SARS-CoV-2 was first detected in 2019 in Wuhan, China, and possibly originated from a recombination event in an ancestor of SARS-CoV-2, a horseshoe bat coronavirus, around 11 years ago via zoonotic transmission from the pangolins [6,7].
Host cell entry is the first step in the viral life cycle. The first step in the life cycle of SARS-CoV-2 is the attachment of the viral particle via the receptor-binding domain (RBD) of its S protein (see below) to the angiotensin-converting enzyme-2 (ACE-2) receptor on the plasma membrane of the pulmonary alveolar epithelial cells and capillary endothelial cells. In some instances, this ultimately leads to severe acute respiratory failure ( Figure 1) [8,9]. Similar to other coronaviruses (CoVs), the size of the SARS-CoV-2 genome is approximately 30 Kb and encodes four main structural proteins, including the Spike glycoprotein (S protein), Envelope (E) protein, Membrane (M) protein, and Nucleocapsid (N protein) [10]. In addition to these four main structural proteins, the SARS-CoV-2 genome codes for sixteen non-structural proteins (NSP 1−16) [11]. Together these proteins facilitate the replication of the virus in the host cell. While effective vaccines and vaccination programs are ongoing, new viral variants are emerging, and infections, hospitalizations and deaths from COVID-19 continue. Consequently, it is of paramount importance to develop effective inhibitors of SARS-CoV-2.

Docking Validation and High Throughput Virtual Screening of BTAs
In this work, a high throughput virtual screening of a library consisting of fourteen marine BTAs, five bromotyrosine compounds with established antiviral properties (1-5, Figure 2), and nine bromotyrosine derivatives (6-14, Figure 3) from the French Polynesian marine sponge, Suberea ianthelliformis, was performed against five SARS-CoV-2 target proteins. The binding potentials of the compounds in this screening library were indicated as S-scores and compared to those of the co-crystallized native ligand for each of the five target proteins for the aim of docking validation and a cross-reference comparison, as shown in Table 1.
For the SARS-CoV-2 Main protease (M Pro ), compound 13 had a better S-score compared to the co-crystallized ligand (S-scores −8.54 vs. −8.25 Kcal/mol, respectively). The other library compounds also demonstrated effective, predicted binding affinities with S-scores between −5.97 and −8.02 Kcal/mol. The S-scores for all fourteen library compounds were greater for the spike glycoprotein (−5.14 to −7.65 Kcal/mol) compared to the spike glycoprotein co-crystallized with a native ligand (S-score = −4.55 Kcal/mol), with compound 3 once more having the best predicted binding affinity (S-score = −7.65 Kcal/mol).
Compound 1 had a similar predicted binding affinity with nucleocapsid phosphoprotein compared to the co-crystallized native ligand (−4.33 vs. −4.44 Kcal/mol, respectively). The other thirteen compounds had higher S-scores (−5.10 to−7.04 Kcal/mol) indicating greater binding affinities, with compound 2 having the highest binding score.   For the SARS-CoV-2 membrane glycoprotein; the spike protein bound to the PD of ACE2 with a dissociation constant of~15 nM and compound 3 having the greatest predicted binding affinity (S-score = −6.28 Kcal/mol).
The 2D interactions of the screening library compounds compared to those of the native co-crystallized ligands of the five SARS-CoV-2 target proteins were studied with binding data for each compound displayed in Table 2.
According to Table 2 Similarly, anomoian F (14) also showed two interactions: a H-bond acceptor between O35 of the ligand and Cys 145 in the protein pocket, and a pi-H interaction between the ligand's six-membered ring and Glu 166. The 2D interactions of the fourteen test compounds from the screening library, in addition to the native ligand with M Pro , are shown in Figure 4.    One interaction: -Pi-H interaction between the ligand's 6-membered ring and Gly 339 in the pocket with distance 3.81 Å and energy score of −0.6 Kcal/mol.
One interaction: -H-bond acceptor between N65 of the ligand and Asn 77 in the pocket with distance 3.56 Å and energy score of −0.9 Kcal/mol.              Fistularin-3 (3) and anomoian C-D (11)(12) showed the best interactions with the membrane glycoprotein as illustrated in Table 2     Compounds 3, 4 and 11 showed the best interactions with the non-structural protein, nsp10, where all formed three interactions with the active pocket, as shown in Table 2 and Figure 8

In Silico Prediction of Pharmacokinetics and Toxicity (ADME/Tox)
The pharmacokinetic properties of the 14 library compounds were calculated in silico using the SWISS-ADME and pkCSM online webtools. These results are shown in Table 3. Twelve of the fourteen compounds had very high logP values (>5) and, consequently, failed to comply with the Lipinski's rule of five requirements in this respect; only compounds 1 and 3 had logP values <5, at 2.44 and 2.97, respectively. Other than 1 and 3, all the other 12 compounds were poorly soluble, and thus had an unfavorable solubility. Compounds 2 and 5 were completely insoluble. Compound 1 had a high solubility and compound 3 was moderately soluble which agreed with their lipophilicity scores. Compound 1, which showed effective pharmacokinetic properties, was the only compound predicted to be able to pass the BBB. Thus, this compound may possibly have side effects in the CNS. Regarding metabolism, all compounds, except 1, 3 and 4, were potential substrates for CYP3A4 enzymes. On the other hand, no compound raised concerns with respect to medicinal chemistry parameters as possibly being pan-assay interference compounds (PAINS). All compounds also showed no potential cardiotoxicity as inhibitors of hERG1.
In conclusion, compound 3 demonstrated the best combination of ADME/Tox properties among all 14 compounds.

Structure-Activity Relationships (SARs)
The SARs of this series of marine alkaloids, based on the results presented in Table 2, are summarized in Figure 9. It seems likely that the presence of the two terminal amines is essential for the interaction with M pro , while the primary terminal amines, in contrast, are not favorable to this interaction. Converting these terminal amines into amides connected to the unsaturated spiro [4,5]decane, as is the case with compound 3, showed the greatest interaction with M pro , having the ability to occupy its four major pockets: S1, S2, S3 and S4. The presence of the tertiary amines is also not favorable, as is the case for compounds 6, 7, 8 and 11. For the interaction, at least one amine must be a secondary amine, as is the case for compounds 9, 10 and 13. However, attaching the amide group to a long saturated aliphatic chain, as in compound 5, provides a better chance of occupying the spike glycoprotein.
Similar to M pro , the presence of the two amides connected to the unsaturated spiro [4,5]decane increases the predicted binding affinities of the nucleocapsid phosphoprotein, membrane glycoprotein, and nsp10, as clearly shown by compounds 3 and 4. However, dissimilar to M pro , the presence of the two terminal amines is not favorable for binding to the nucleocapsid phosphoprotein, membrane glycoprotein, or nsp10, as is the case for compounds 11 and 12. The presence of a terminal hydroxyl group on the other side of the compound increases the binding interactions in these compounds. The presence of alpha-beta unsaturated compounds, such as those in 6-10, decreases binding ability. As compound 3 was the only compound that showed a combination between effective ADME/Tox properties and interactions with high S-scores in all five SARS-CoV-2 target proteins, it was selected for conducting a 100 ns MD simulation with the five target proteins. Figures 10 and 11 show the RMSD fluctuations of protein-ligand complexes with respect to the initial structure, and the radius of gyration, respectively, for compound 3 with the five targets, respectively. This enabled the analysis of the stability of the simulated system throughout the 100 ns MD simulations. As expected, all complexes showed predicted, small RMSD fluctuations within only 2 Å, confirming their high stability throughout the whole simulation, where compound 3 showed the greatest stability with membrane glycoprotein (PDB ID: 6M17). Moreover, the radius of gyration was also consistent with high stability, with all fluctuations being within 0.05 nm.

Preparation of Protein Structures
The X-ray crystal structures for the five target proteins from SARS-CoV-2, including the main protease (M Pro ; PDB ID: 6LU7), spike glycoprotein (PDB ID: 6VYB), nucleocapsid phosphoprotein (PDB ID: 6VYO), membrane glycoprotein (PDB ID: 6M17), and nonstructural protein 10 (nsp10;PDB ID: 6W4H), were retrieved from the Protein Data Bank (http://www.pdb.org, accessed on 1 July 2021). Their resolutions were 2.16 Å, 3.20 Å, 1.70 Å, 2.90 Å, and 1.80 Å, respectively. All water molecules were removed from these crystal structures with only main-chain amino acids retained. An AMBER (AMBER10:EHT) force field was used for the energy minimization of these five X-ray crystal structures using parameters suitable for proteins and nucleic acids (ff10) and small molecules (EHT). Protons were added by employing the 3D protonation feature in MOE v.2019.01; Asn, Gln and His flips were allowed during 3D protonation. Complexes were then refined to a RMS gradient of 0.1 Kcal/mol/Å.

In Silico Prediction of Pharmacokinetics and Toxicity
The pharmacokinetic properties of the fourteen compounds in our screening library were calculated using the SWISS-ADME webtool (https://www.swissadme.ch, accessed on 28 June 2021). The properties predicted here were lipophilicity, reported as Log Po/w (WLOGP); water solubility class; and Blood-brain barrier (BBB) penetration, in addition to medicinal chemistry parameters employing pan-assay interference alerts (PAINS) [39,40]. Additionally, the potential toxicity profiles of these molecules were predicted using the pkCSM online webtool (http://biosig.unimelb.edu.au/pkcsm/prediction, accessed on 28 June 2021) to predict the safety of these small molecules upon ingestion in human and animal models, with respect to toxicological effects on hERG-I inhibition [41].

Molecular Dynamics Simulation for Compound 3
Compound 3, Fistularin-3, displayed the best binding interactions and free energies with the five SARS-CoV-2 target proteins among the fourteen library compounds investigated. It also showed the best pharmacokinetics properties. Accordingly, it was subjected to 100 ns molecular dynamics investigation against the five SARS-CoV-2 target proteins. MD simulations were performed using the GROMACS 2021 software package with the CHARMM36 force field used for protein topology preparation and the official CHARMM General Force Field server (CGenFF) used for ligand topology preparation. The solvation method used was a dodecahedron box of common simple point charge (SPC) water model with explicit solvent periodic boundary conditions. Charge neutralization using sodium and chloride ions was performed for the five solvated complexes. These systems were subjected to energy minimization to resolve steric clashes or inappropriate geometry employing the steepest descent method of 5000 steps. System equilibration was also set to ensure a reasonable starting structure using NVT and equilibration under constant number of particles, volume, and temperature (NVT) for 100 ps using a Berendsen thermostat [42]. Then, re-equilibration was performed for another 100 ps under constant pressure (Isothermal-isobaric (NPT) ensemble) using the Parrinello-Rahman barostat using a time step of 2 fs for each equilibration round [43]. Finally, an MD production phase was performed for 100 ns using a time step of 2 fs at a constant temperature of 300 K and constant pressure of 1 atm. Simulation results were analyzed using Visual Molecular Dynamics (VMD) software, ver.1.9.3 [44].

Post MD Analysis, Trajectory Post-Processing and MM/PBSA Calculations
After determining the trajectories of the five complexes resulting from the MD simulation of compound 3, the complexes were re-centered and rewrapped within unit cells using the trjconv function of GROMACS. The stabilities of trajectories were then determined throughout the 100 ns simulation using the radius of gyration and the root-mean-square deviation (RMSD) of the protein backbone referenced to its initial position at 10 ps intervals. Lastly, g_mmpbsa was employed using Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) binding free energy [45] to calculate relative binding free energies according to the following equation: ∆G bind = ∆E gas + ∆G solvation − T∆S ∆E gas = E int + E vdw + E elec (3) G solvation, GB = G GB + G nonpolar, solvation − G ligand (5) ∆G nonpolar = γSASA + β (6)
Supplementary Materials: The following are available online. Table S1: Final docking validation, Figure S1: 3D interactions between compounds 1-14 and the native ligand with MPro (PDB ID: 6lu7), Figure