Abrogation of SARS-CoV-2 interaction with host (NRP1) neuropilin-1 receptor through high-affinity marine natural compounds to curtail the infectivity: A structural-dynamics data.

The evolution of new severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) variants around the globe has made the coronavirus disease 2019 (COVID-19) pandemic more worrisome, pressuring the health care system and resulting in an increased mortality rate. Recent studies recognized neuropilin-1 (NRP1) as a key facilitator in the invasion of the new SARS-CoV-2 into the host cell. Therefore, it is considered an imperative drug target for the treatment of COVID-19. Hence, a thorough analysis was needed to understand the impact and to guide new therapeutics development. In this study, we used structural and biomolecular simulation techniques to identify novel marine natural products which could block this receptor and stop the virus entry. We discovered that the binding affinity of CMNPD10175, CMNPD10017, CMNPD10114, CMNPD10115, CMNPD10020. CMNPD10018, CMNPD10153, CMNPD10149 CMNPD10464 and CMNPD10019 were substantial during the virtual screening (VS). We further explored these compounds by analyzing their absorption, distribution, metabolism, and excretion and toxicity (ADMET) properties and structural-dynamics features. Free energy calculations further established that all the compounds exhibit stronger binding energy for NRP1. Consequently, we hypothesized that these compounds might be the best lead candidates for therapeutic interventions hindering virus binding to the host cell. This study provides a strong impetus to develop novel drugs against the SARS-CoV-2 by targeting NRP1.

Introduction

The global coronavirus disease 2019 (COVID-19) pandemic was caused by acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a highly contagious virus first identified in Wuhan, China [1]. In 2013, a similar coronavirus named SARS-CoV caused an epidemic. SARS-CoV-2 infects the lower respiratory system and spreads quickly, mainly via pharyngeal viral shedding [2,3]. Angiotensin-converting enzyme 2 (ACE2), a host receptor, is responsible for the uptake of SARS-CoV and SARS-CoV-2 viruses and is thus considered a primary target for viral entry [4,5].

Recently, neuropilin-1 (NRP1) was identified as an additional entry channel that facilitates the entry of SARS-CoV-2 into the host cell. NRP1 is a member of the catalytic and signalling proteins family that is reported to be involved in the cellular invasion of SARS-CoV-2 and to potentially cause in vitro infectivity of SARS-CoV-2 [6]. However, it is uncertain whether NRP1 can enable the receptor-mediated endocytosis of the virus. NRP1 is also recognized as a cellular signalling and cell surface receptor [7,8]. The host receptor binding protein of SARS-CoV-2, known as the spike protein, is cleaved into S1 and S2 polypeptides by protease and furin within the host cell, thus exposing the CendR motif in the cleaved S1 protein. The CendR motif is named after the C-end terminal rule requiring the presence of a cationic amino acid, generally arginine, at the carboxyl terminus, which causes an RXXR configuration in the ligand. The binding pocket of CendR is present in NRP1's b1 domain [9]. Daly et al. found that the virus's infectivity is caused by the binding of the CendR motif to NRP1 in the S1 protein. Neutralization of the virus is initiated by the attachment of a monoclonal antibody to the extracellular binding pocket (ECP) of the CendR region of the b1 domain in NRP1 [10]. The CendR motif b1 domain is exposed as soon as the viral spike protein is cleaved by the proteases and thus mediates the binding and entry through NRP1 [11]. This motif's C-end terminal rule name is due to the interspersed cationic amino acids, particularly arginine, that result in an RXXR configuration [[12], [13], [14]]. According to Daly et al. [7], the CendR motif in the SARS-CoV-2 S1 protein interacts with NRP1 and increases viral infectivity. NRP1 RNA expression analysis of bronchial and alveolar cells isolated from patients with severe COVID-19 revealed elevated expression of NRP1 in SARS-CoV-2-positive cells but not in uninfected cells [8]. Blocking the association between NRP1 and SARS-CoV-2 may be a useful therapy in the battle against COVID-19, but further research is required. The results of these studies show that in addition to the function of ACE2 in causing the cellular entry of SARS-CoV-2, NRP1 can also perform as a mediator to increase the virus's infectivity.

The conclusions of these well-designed studies indicate that NRP1 may act as a host cell mediator to enhance the virus's attachment and entry into host cells; thus, it may maximize tissue tropism of SARS-CoV-2 [15]. Therefore, in the present study, molecular modelling-based approaches were employed to target NRP1 in order to identify novel hits against SARS-CoV-2 within a comprehensive marine natural products (CMNP) database (https://www.cmnpd.org/). Computational structural biology and biophysical approaches are the most widely used approaches in drug discovery for the treatment of viral pathogens [[16], [17], [18], [19], [20], [21], [22], [23], [24], [25]]. We used virtual drug screening approaches and binding free energy calculations to decipher the binding, possible mechanism, and long-term association strength between the identified drugs and NRP1. These efforts to discover effective drug candidates for the treatment of COVID-19 will help combat the infection caused by SARS-CoV-2.

Materials and methods

Marine library filtration

Before structure-based VS, the CMNP database was filtered based on drug-like soft, toxicophores, removal of Pan Assay Interference (PAINS) and Eli Lilly MedChem rules in FAFDrugs4 server (https://fafdrugs4.rpbs.univ-paris-diderot.fr/) to obtain drug-like, non-toxic, and non-PAINS compounds [26]. The different parameters used in discarding non-drug-like molecules are presented in Table 1 .

Table 1

Different drug-like parameters used in CMNP filtration.

Parameters	Value
Molecular weight	100–600
logP	−3 to 6
Hydrogen bond acceptor	≤12
Hydrogen bond donors	≤7
tPSA	≤180
Rotatable Bonds	≤11
Rigid Bonds	≤3
Rings	≤6
Max size system ring	≤18
Carbons	3–35
HeteroAtoms	1–15
H/C Ratio	0.1 to 1.1
Charges	≤4
TotalCharge	−4 to 4
RO5 Violations	No

Protein and ligand structure preparation

The crystal structure of NRP1 (PDB ID: 1KEX) was retrieved from Protein Databank (http://www.rcsb.org/) [27]. The structure was prepared by removing water molecules and co-crystallized ligands using the protein preparation wizard incorporated in Schrödinger software [28]. The Amber force field was used for protein's energy minimization using Chimera [28,29]. For ligands preparation, the filtered compounds were imported to PyRx and subjected to energy minimization using MM2 force field [30].

High-throughput VS

AutoDock Vina was used to virtually screen compounds at the binding site of protein collectively formed by residues Tyr297, Trp301, Thr316, Asp320, Ser346, Thr349 and Tyr353 [31]. The grid box size was defined as 6.04 x −70.97 × 24.028, while the grid dimensions of 56 Å × 58 Å × 50 Å were generated. The aforementioned experimentally validated residues were selected to specify the final active site grid and the maps generated by the AutoDock tool. Three steps of VS were performed to screen the CMNP database [32]. Initially, the whole database was screened using AutoDock Vina with exhaustiveness set as 16, and then the top hits were screened again using exhaustiveness as 32. Finally, for the best scoring compounds, induced-fit docking(IFD) was carried out using 64 exhaustiveness to confirm the final hits. For IFD, we used AutoDockFR–AutoDock for Flexible Receptors (ADFR) software which handles the receptor sidechain conformational optimization up to 14 different sidechains to enhance the success rate of docking [33]. AutoDockFR achieved higher accuracy than AutoDock Vina in cross-validation, and also the speed of the docking is much higher.

Absorption, distribution, metabolism, excretion and toxicity (ADMET) properties and bioactivity prediction

For the top hits, parameters like drug-like, lead likeness, toxicity, pharmacokinetics and solubility were also explored using SwissADME [34]. Molinspiration, a cheminformatic tool (https://www.molinspiration.com/cgi-bin/properties), was used to foresee the IC50 value of each compound. For instance, Molinspiration is widely employed by thousands of studies to forecast bioactivity scores.

Determination of dissociation constant (KD.)

Furthermore, to give a persuasible understanding of the dissociation constant (KD), PROtein binDIng enerGY prediction (PRODIGY) (https://wenmr.science.uu.nl/prodigy/) server was used to computationally predict KD. values for different biological complexes [35].

Molecular dynamics simulation of protein-ligand complexes

All-atoms MD simulation of the top hits from the CMNP database was performed using the AMBER18 package [36]. For drug topologies, antechamber module was used while the Amber general force field (GAFF) and ff14SB forcefields were employed for parameterization of the ligands and NRP1, respectively. To solvate each protein-ligand complex, a TIP3P water box was used, while Na + counter ions were used to neutralize each system subsequently. Energy minimization of systems was carried out in two stages (steepest descent and conjugate gradient), followed by heating and equilibration [23]. The Particle Mesh Ewald (PME) algorithm was used to treat long-range electrostatic interactions [37]. A 1.4 nm cutoff value was set for Van Der Waal interactions, as well as for Columbic interactions of short-range. Langevin thermostat was employed to temperature constant at 300 K, whereas for pressure control, Berendsen barostat was considered. A time step of 2fs and a total simulation time of 200 ns for each complex was performed. The dynamics, stability and other features of the ligand-protein complexes were evaluated by using CPPTRAJ and PTRAJ [38].

The binding free energy calculations

For all protein-ligand complexes, the free binding energy was calculated using the script MMPBSA.PY by considering 500 snapshots using equations described below [[39], [40], [41], [42]]. Different studies have used this free energy calculation method to estimate the binding of ligands for a given biological macromolecule [43,44].

Here, ΔGbind denotes total free binding energy, while others denote the free energy of the protein, the ligand and complex. The following equation was used to calculate specific energy term contribution to the total free energy:

Bonded, electrostatic, polar, non-polar and van der Waal energy terms are represented by the above equation.

Results and discussion

Computer VS and binding modes of top hits

The outbreak of SARS-CoV-2 has distressed the world by spreading exponentially to every nook and cranny of the world. The pathophysiology of this virus shows that it primarily affects the respiratory tract by entering the host cell through Spike protein. The attachment of the Spike receptor-binding domain to the host ACE2 protein enables the entry to the host cell, and it is considered the main drug target for novel drug discovery [45,46]. However, recent findings deemed NRP1 as a supplementary cellular intermediary that enables SARS-CoV-2 attachment and entrance into host cells [7,8,11]. Therefore, it could serve as an alternative therapeutic target. Keeping in view the importance of NRP1 in the entry of SARS-CoV-2 and infection herein, we applied a computational VS approach to identify potential inhibitors from the marine CMNP database using AutoDock Vina. The database has 47,451 molecules collected from different marine sources. Our analysis identified that among the 47,451 compounds, only 21,000 compounds passed the ADMET criteria. Using AutoDock Vina for virtual screening reported that for these 21,000 compounds, the scores ranged from −8.0 to −2.3 kcal/mol. Compounds with a score lower than −6.0 kcal/mol were selected for further analysis. This was done to select highly efficient binders as a very small number of compounds from the entire library were found to have stable binding energy with the protein because of weak interaction profile. Therefore, to increase selection chances of efficient binding molecules to the protein, a stringent criteria of −6.0 kcal/mol to −8.0 kcal/mol was used. With this criterion, only those compounds were selected which interact with the key hot spot residues of the receptor protein. Based on the high docking score and interactions with the key residues, ten compounds were shortlisted for further analysis. The identified top ten compounds are given in Table 2 . These compounds were further analyzed for drug-like rules, which include Lipinski, Ghose, Veber, Egan and Muegge parameters for drug feature calculation [47]. The shortlisted compounds include CMNPD10175, CMNPD10017, CMNPD10114, CMNPD10115, CMNPD10020. CMNPD10018, CMNPD10153, CMNPD10149 CMNPD10464 and CMNPD10019. The shortlisted compounds were reported to act drug-like with no toxicity and good pharmacokinetic properties. From the top ten compounds, only five compounds, presented in Fig. 1 , were selected for interactions, molecular dynamics simulation and free energy calculations. The following sections describe the binding mode of each of these compounds in more detail.

Table 2

Top ten marine drug compounds shortlisted through VS and applied drug-like filters. The docking score of each compound is given in kcal/mol.

Compound ID	Compound Names	drug like	Lead likeness	Toxicity	Pharmacokinetics	Solubility	Score
CMNPD10175	Madangolide	Yes	Yes	No	Good	Good	−8.0
CMNPD10017	Isogranulatimide	Yes	Yes	No	Good	Good	−6.9
CMNPD10114	Luisol B	Yes	Yes	No	Good	Good	−6.9
CMNPD10115	Luisol A	Yes	Yes	No	Good	Good	−6.8
CMNPD10020	Meridianin A	Yes	Yes	No	Good	Good	−6.6

Fig. 1

2D structure representation of the top five hits along with their CMNPD. database accession ID and original scientific names are given. The structural scaffold of all the compounds seems very similar, which reflects the same activity of these compounds.

Binding mode of CMNPD10175

Madangolide (CMNPD10175) is a marine natural compound originally isolated from cyanobacteria known Lyngbya bouillonii, which entails one large macrocyclic lactone ring [48]. It belongs to a group of antibiotics known as macrolides. Herein, Madangolide has been identified as the most potent drug with a docking score of −8.0 kcal/mol. This compound formed two hydrogens and two pie-alkyl bonds with the key active site residues. Among the two hydrogen bonds, one was formed with Phe297 while the other one was reported with Lys351. The residue Phe297 was also involved in pie-alkyl along with Tyr353. The interaction pattern of Madangolide (CMNPD10175) is given in Fig. 2 a.

Fig. 2

Binding modes of Madangolide (CMNPD10175) (a) and Isogranulatimide (CMNPD10017) (b). The hydrogen-bonding interactions are shown in green colour, while the pie-alkyl interactions are shown as pink.

Binding mode of CMNPD10017

Isogranulatimide (CMNPD10017), marine pyrrolocarbazole, is a cell-permeable alkaloid primarily used for phosphorylation and dephosphorylation and has been reported to act as competitive inhibitors in cancer [49]. Initially isolated from a marine source, Didemnum conchyliatum, this compound act as a major G2 cell cycle checkpoint inhibitor. This compound identified herein through virtual drug screening possesses a docking score of −6.9 kcal/mol against the NRP1. Similar to the Madangolide, this compound formed 4 key interactions with the NRP1 receptor. Among the total interactions two hydrogen bonds with Asp320 and Thr413 while the two pie-alkyl interactions were formed with Tyr297, which is considered as an important residue for the recognition. The interaction pattern of Isogranulatimide (CMNPD10017) is given in Fig. 2b.

Binding mode of CMNPD10115

Likewise, the CHMPDB10115, Luisol B (CMNPD10115) also belongs to Naphthofurans class of marine natural products produced by a marine bacterium from the Streptomyces (Actinomycetales) genus [50,51]. This compound identified herein through virtual drug screening possesses a docking score of −6.9 kcal/mol against the NRP1. Unlike the Madangolide and Isogranulatimide, this molecule formed three hydrogen bonds and two pie-alkyl interactions with the key residues. Among the hydrogen bonds Tyr297, Tyr353 and Thr413 are involved, while the two pie-alkyl interactions are formed with Tyr297 and Trp301. This compound shows robust interactions with the key active site residues. The interaction pattern of Luisol B (CMNPD10114) is given in Fig. 3 a.

Fig. 3

Binding modes of Luisol B (CMNPD10114) (a) and Luisol A (CMNPD10115) (b). The hydrogen bonding interactions are shown in green colour while the pie-alkyl interactions are shown as pink.

Binding mode of CMNPD10114

Luisol A (CMNPD10114) belongs to Naphthofurans class of marine natural products produced by a marine bacterium from the Streptomyces (Actinomycetales) genus [50,51]. This compound secured a docking score of −6.8 kcal/mol against the NRP1. This molecule formed two hydrogen bonds and two pie-alkyl interactions with the key residues. Residues Trp301 and Glu348 are involved in hydrogen bonding while the two pie-alkyl interactions are formed with Tyr297 and Trp301. The interaction pattern of Luisol A (CMNPD10115) is given in Fig. 3b.

Binding mode of CMNPD10020

Meridianin A (CMNPD10020) is a marine-derived hydroxyindoles primarily identified as kinase inhibitor [52]. Meridianin A identified as a potential drug with a docking score −6.6 kcal/mol. This compound formed two hydrogens and five pie-alkyl bonds with the key active site residues. Among the two hydrogen bonds were formed with Trp413 while the five pie-alkyl interactions involved Tyr297, Trp301 and Tyr353 residues. The interaction pattern of Meridianin A (CMNPD10020) is given in Fig. 4 .

Fig. 4

Binding modes of Meridianin A (CMNPD10020). The hydrogen-bonding interactions are shown in green colour, while the pie-alkyl interactions are shown as pink.

For instance, a comparative analysis of these top five hits with the already experimentally reported compounds revealed that the shortlisted compounds possess better docking scores than the experimentally reported. The experimentally reported best compound EG00229 has docking score of −6.4 kcal/mol, while the others have docking energy as follow: EG01377 (docking score = −5.8 kcal/mol) [53], CB7739526 (docking score = −4.1 kcal/mol) [54] and NRPa308 (docking score = −3.5 kcal/mol) [55]. In contrast, the docking scores of our shortlisted compounds are far better than the experimentally reported, which confirms the stronger inhibitory effects of our compounds. Similarly, a recent study [56] also reported some compounds including Amantadine, Carvacrol, Thymoquinone, and Thymol identified through VS against the NRP1 also possess lower docking scores than our reported top hits. This shows the robustness of our top hits, which bind more strongly than the previously reported.

Prediction of bioactivity and KD. Determination

On the other hand, molinspiration predicted the bioactivity of each compound against each class of receptors. The server's results demonstrate that all of the compounds on the shortlist are active against the nuclear receptor target ( Table 3 ). The reported scores for CMNPD10175 (0.34), CMNPD10017 (0.64), CMNPD10114 (0.50), CMNPD10115 (0.46) while CMNPD10020 possess (0.04) bioactivity score against the nuclear receptor target. The scores for other classes such as GPCR, Ion Channels, Kinases, protease and enzymes are also given in Table 3. As a result, these findings conclude that the compounds shortlisted could effectively inhibit NRP1 in the experimental design and could be evaluated in clinical trials.

Table 3

Top five marine drug compounds shortlisted through VS and drug-like rules. The docking score of each compound is given in kcal/mol.

Compound ID	GPCR ligand	Ion Channels	Kinase Inhibitors	Nuclear Receptor ligand	Protease inhibitor	Enzymes inhibitor	Prodigy score
CMNPD10175	0.18	−0.06	−0.42	0.34	0.28	0.21	−7.9
CMNPD10017	−0.07	−0.10	1.29	0.64	−0.20	0.41	−7.3
CMNPD10114	0.25	0.17	−0.10	0.50	0.21	0.51	−7.3
CMNPD10115	−0.04	0.01	−0.20	0.46	−0.11	0.34	−7.3
CMNPD10020	0.20	0.34	1.0	0.04	−0.46	0.60	−7.1

We also determined the KD of all the complexes to provide a clear insight into the binding variations. KD is a method for assessing and ranking the strength of biomolecular interactions [57]. This method is widely used for the binding strength determination of biological molecules [58,59]. We employed PRODIGY, to calculate the binding affinity KD for these protein-drugs complexes. Prodigy reported the binding score for each drug as −7.9 for CMNPD10175, -7.3 for CMNPD1017, CMNPD10114 and CMNPD10115, respectively. While for CMNPD10020, the prodigy score was predicted to be −7.1 kcal/mol. This shows the tighter binding of these molecules with the host receptor NRP1 and its inhibition. The predicted scores by prodigy are given in Table 3.

Deciphering structural stability via RMSD

Predictions made by molecular docking about enzyme-ligand docked conformation stability were tested in a dynamic environment over 200 ns of MD simulation time. To understand the generated MD simulation trajectories, CPPTRAJ. module of AMBER was employed to carry out several statistical parameters that helped to understand overall complexes structural stability. In this regard, the first root mean square deviation (RMSD) was performed that plotted all Carbon alpha atoms deviation from original docked conformation versus time. Higher RMSD. implies higher structural fluctuations with respect to a reference, whereas lower RMSD. corresponds to good structure stability, i.e. higher intermolecular affinity and better-docked mode. The mean RMSD of CMNPD10017, CMNPD10020, CMNPD10114, CMNPD10115, and CMNPD10149 was reported to be 1.2 Å, 1.4 Å, 1.3 Å, 1.3 Å and 1.2 Å, respectively. Overall, all the complexes indicated good structural stability as no significant spike in the RMSD was highlighted. The CMNPD10149-enzyme complex remained quite stable till 50 ns with RMSD. ~1 Å. The followed RMSD suffers from an initial minor deviation touching RMSD of 2 Å. Then, the RMSD plot is more uniform, with no global or local structural changes noticed until 170 ns. Afterwards, a sudden rise and fall of RMSD till 175 ns was depicted. Towards the end, the RMSD of this system was seen as the highly stable indicating equilibrium of the intermolecular interactions and overall system. The CMNPD10017-enzyme complex is subject to a continuous pattern of RMSD with a minor steady increase till 165 ns. The highest RMSD seen for this system is 2.5 Å. Towards the end time, system RMSD is drop-down touched 1.5 Å and suffered from a sudden upsurge. This concluded the system RMSD not converged and need more MD simulation time. The CMNPD10114 system is more consistent with no structural alterations and a mean RMSD of 1.3 Å to 1.4 Å. CMNPD10115 and CMNPD10020 complexes like CMNPD10114 behave very stably except for one sharp spike at 125 ns in the case of CMNPD10115. The last three systems were reported to show the greater intermolecular affinity of the compounds for the enzyme that further suggested that proposed docked complexes by docking studies are valid in terms of lowest energy minima and conformational stability. The RMSDs of all the complexes are given in Fig. 5 .

Fig. 5

Thermodynamics stability calculated as RMSD of all the complexes. (a) CMNPD10175, (b) CMNPD1017, (c) CMNPD10114, (d) CMNPD10115 and (e) CMNPD10020 represent each complex.

Investigating residues fluctuations

Investigating the residues fluctuation, we employed root-mean-square fluctuation (RMSF) to calculate the residual flexibility of carbon alpha residues of the targeted systems. Generally speaking, all the systems mean RMSF is very low, indicating the residues of enzyme enjoying excellent energy minima in the presence of compounds. This can be indeed affirmed by the very low RMSF of enzyme binding site residues which is well below under 2 Å. However, the C- and N-terminal of the enzyme were found highly unstable, which might be flexible and is a mechanism to attain a more stable conformation of the compounds at the active pocket. RMSF data complement RMSD as both agree on the highly stable nature of complexes. The RMSFs of all the complexes are given in Fig. 6 .

Fig. 6

Residual flexibility calculated as RMSF of all the complexes. All the complexes CMNPD10175, CMNPD1017, CMNPD10114, CMNPD10115 and CMNPD10020 represented with different colour.

Determining strength of intermolecular interactions

In biological systems, hydrogen bonding is key to strong intermolecular binding and therefore easing molecular recognition and ultimately performing biological functions. All complexes were subjected to hydrogen bond analysis to calculate the number of hydrogen bonds formed in each frame of MD simulation. As can be visualized in the figure, all the targeted complexes are dominated by a rich network of hydrogen bonds reflecting on the strong stability of the docked compounds with the enzyme. The average number of hydrogen bonds in each frame of the complexes is 85, 83, 83, 82 and 82, respectively. As observed in the above RMSD and RMSF, the systems are showing significant stability, which is the outcome of a large number of hydrogen bonding. The H-bonds count of all the complexes is given in Fig. 7 .

Fig. 7

Thermodynamics stability calculated as RMSD of all the complexes. (a) CMNPD10175, (b) CMNPD1017, (c) CMNPD10114, (d) CMNPD10115 and (e) CMNPD10020 represent each complex.

Unveiling atomic level interaction energies

Determining binding free energies at the residue level of protein-ligand complex is now in the routine of computer-aided drug designing (CAAD) to cross-validate docking and MD simulation results. This was achieved by performing MM-GBSA analysis on the simulation trajectories by picking all the frames. Binding free energies estimation by MM-GBSA is a common practice as they are more reliable than traditional docking scores and less computationally expensive compare to Alchemical free energy methods. As can be understood from the data tabulated in Table 4 , the CMNPD10017 is relatively more stable by securing the lowest binding energy of −43.08 kcal/mol, followed by CMNPD10020 (−40.86 kcal/mol). This can be reasoned by the structural deviations the complexes acquired during simulation time; these deviations might help these two complexes to get more stable conformation by establishing new chemical interactions. The last three complexes, CMNPD10114, CMNPD10115 and CMNPD10149, are classified as less stable than the mentioned before that might be because of the relaxed complex structure. The complexes preferred to deviate from the original conformation as this might lose their significant bonding and will push the systems to instability. CMNPD10017, CMNPD10020 and CMNPD10149 interactions with the enzyme are dominated by van der Waals energy as well as a favourable contribution from electrostatic energy was also unraveled. In the case of CMNPD10114 and CMNPD10115, the systems favoured electrostatic energy to produced stable complex while less significant help was considered from van der Waals energy. In all complexes, the non-polar energy of the solvation phase seems also to contribute positively, making the systems stable. The polar energy is a non-contributor to system stabilization as can be interpreted by positive values. Together, docking, MD simulations and binding free energies all are in support of good docked conformation of the compounds to the enzyme and resulting in highly stable complexes.

Table 4

Binding free energy results of all the complexes calculated using the MM-GBSA approach. The free score of each compound is given against it in kcal/mol.

Complex	vdW	Electrostatic	GB	ESURF	Total
CMNPD10017	−30.97	−25.60	16.78	−3.29	−43.08
CMNPD10020	−27.00	−18.17	7.08	−2.59	−40.86
CMNPD10114	−24.65	−93.53	89.10	−3.30	−32.38
CMNPD10115	−20.76	−140.92	131.88	−2.99	−32.80
CMNPD10149	−27.93	−11.65	4.83	−3.00	−34.45

Conclusions

In conclusion, this study identified novel hits from the marine sources against SARS-CoV-2 NRP1 receptor. Though the findings of the current work are promising, the study suffers from several limitations of the methods applied herein. Docking limitations are overcome by the long run of MD simulations and downward sophisticated computational analysis. We further determined the bioactivity and KD of the top hits. CMNPD10175, CMNPD10017, CMNPD10114, CMNPD10115, CMNPD10020. CMNPD10018, CMNPD10153, CMNPD10149 CMNPD10464 and CMNPD10019 which were identified through VS and further explored by analyzing its ADMET properties and structural-dynamics features. Free energy calculations further established that all the compounds exhibit strong binding energy for NRP1. Though the findings of the current work are promising, the study suffers from several limitations of the methods applied herein. For example, the docking predictions are often found non-reliable. The force fields used in MD simulations need further refinement whereas the binding free energy methods lacked entropy estimation as well as missed the role of water molecules that bridge ligands with receptor hotspot residues. Consequently, we hypothesized that these compounds might be the best lead candidates for therapeutic interventions. This study provides a strong impetus to develop novel drugs against the SARS-CoV-2 by targeting NRP1.

Author contributions

“Conceptualization, AK and A.M.; methodology, A. K; software, DW; validation, DW, FA and JA; formal analysis, AK DW SA; investigation, KK, SS.A. JA, NU, AS; resources, DW; data curation, AK, AS; writing—original draft preparation, AK, KK, SS.A. JA, NU, AS; writing—review and editing, A.M, DW; visualization, AKFA, JA; supervision, DW; project administration, DW; funding acquisition, DW All authors have read and agreed to the published version of the manuscript.”

Funding

Dong-Qing Wei received funding from Key Research Area Grant 2016YFA0501703 of the “Ministry of Science and Technology of China”, “the ”. The grants Numbers are 32070662, 61832019, 32030063. This work is also supported by “” under the grant number 19430750600 and Joint Research Fund and Joint Research Funds for Medical and Engineering and Scientific Research from grant number YG2021ZD02.

Data availability statement

All the data is available on RCSB, UniProt and any simulation data would be provided on demand. The accession numbers to access this data are given in the manuscript.

Declaration of competing interest

Declared None.