In-silico natural product database mining for novel neuropilin-1 inhibitors: molecular docking, molecular dynamics and binding energy computations

Abstract

In the search for new metabolite inhibitors, a natural product activity and species source (NPASS) database was virtually screened using AutoDock software to identify potential NRP1 inhibitors. NPASS compounds complexed with NRP1 were subjected to molecular dynamics (MD) simulations. Furthermore, NPASS-NRP1 binding affinities were calculated using the MM-GBSA approach. Based on calculated binding energies, kamolonol (NPC146388) demonstrated lower NRP1 binding affinity than the co-crystallized HRG/Arg-1 ligand with binding energy (ΔG_binding) values of –34.5 and –32.0 kcal/mol, respectively. Structural and energetic analysis showed high stability for kamolonol and HRG/Arg-1 with NRP1 over the 200 ns MD simulations. The studied physicochemical properties of kamolonol and HRG/Arg-1 revealed that these compounds obey Lipinski's rule of five. ADMET characteristics of kamolonol and HRG/Arg-1 were predicted, and kamolonol showed better ADMET properties compared to HRG/Arg-1. Based on these results, kamolonol is a promising NRP1 inhibitor worthy of further experimental assays as anti-carcinoma remediation.

KEYWORDS:

• Cancer
• neuropilin-1 (NRP1)
• docking calculations
• MD simulations
• pharmacokinetic study

1. Introduction

Biological alteration that results in unregulated cell growth [1] is referred to as cancer and is a worldwide health challenge with poor clinical outcomes [2,3]. Approximately 10 million deaths were reported in 2020 by global cancer statistics [4]. Despite ongoing progress in surgery and radiotherapy, patients with advanced tumours still have poor prognoses [5]. Moreover, there is only a small number of newly approved cancer medications released annually by the US Food and Drug Administration (FDA) [6]. With resistance to therapeutic drugs a substantial problem in treating cancer [7], there is a critical need for novel cancer therapies.

Numerous cancers, including gastric, urinary, bladder, breast, lung, and esophageal, have high levels of neuropilin-1 (NRP1) protein [8–11]. NRP1 is a cell surface receptor that binds vascular endothelial growth factor (VEGF), a potent mediator of endothelial permeability, chemotaxis, and proliferation. It is also a regulator for tumour development, proliferation, invasion, metastases, and prognosis [12–16]. NRP1 enhances angiogenesis by forming a direct complex with endothelial growth factors A and 2 (VEGFA and VEGFR2) [17]. It also promotes the tumour's development after binding with the VEGFA due to its ability to promote RhoA activation [18]. Besides, NRP1 has the ability to accelerate neoplasm progression by steadying the role of regulatory T cells (Tregs) and banning the tumour-associated macrophages (TAM) from penetrating normoxic cancer areas [19]. Although NRP1 protein has received considerable critical attention as a potential therapeutic target for cancer therapy [20,21], no drug has been approved as an NRP1 inhibitor.

Natural products (NPs) retrieved from microbes and plants have a vital role in drug discovery, particularly against carcinoma and contagious illnesses [22,23]. Several drugs are natural products from plant sources, such as paclitaxel, morphine, and vinblastine [24]. Therefore, there is abundant opportunity for mining novel anticancer drug candidates as NRP1 inhibitors from the Natural Product Activity and Species Source (NPASS) database. In the current study, NPASS database was screened for promising NRP1 inhibitors utilizing molecular docking computations [25]. Subsequently, promising NPASS compounds were selected and subjected to molecular dynamics (MD) simulations and the corresponding binding energies were estimated using the MM-GBSA approach. The stability and tightness of the most potent NPASS compounds within NRP1 were selected. Thus, the current study highlights the potentiality of NPASS candidates as NRP1 inhibitors and illustrates a promising drug candidate for clinical consideration.

2. Computational methodology

2.1. NRP1 preparation

The neuropilin-1 b1 domain (NRP1) atomic coordinates in complex with L-homoarginine (HRG/Arg-1) were retrieved from the protein databank (PDB code: 5ijr; resolution: 1.52 Å [26]) and considered as a template for all in-silico computations. The crystallographic water molecules, ions, and ligand were eliminated from the NRP1 protein, keeping only the amino acids. The protonation state of the amino acids was investigated using the H++ server, and all missing hydrogen atoms were inserted [27].

2.2. Database preparation

The Natural Product Activity and Species Source (NPASS) database, containing more than 35,000 natural compounds, was downloaded in SDF format [25]. The duplicated compounds were eliminated based on the International Chemical Identifier (InChIKey) [28]. The three-dimensional (3D) structures of the NPASS compounds were generated using Omega2 software [29,30] and minimized utilizing the MMFF94S force field integrated inside SZYBKI software [31,32]. The protonation state of the NPASS compounds was examined using fixpka algorithm within the QUACPAC software, respectively [33]. The Gasteiger-Marsili method was used for the determination of the partial charges of NPASS compounds [34]. The prepared files are present in a CompChem database (www.compchem.net/ccdb) that is publically accessible. The schematic diagram of the utilized computational approaches is depicted in Figure 1.

Figure 1. The schematic diagram of the utilized computational approaches in the virtual screening process of the Natural Product Activity and Species Source (NPASS) database.

2.3. Molecular docking

AutoDock Vina1.1.2 software was used for executing molecular docking calculations [35]. The NRP1 protein was converted into pdbqt format using MGL tools (version 1.5.7) [36]. The details of the employed protocol for the docking calculations are described elsewhere [37–39]. Except the exhaustiveness parameter, which was set to 50 and 200 in fast and expensive calculations, respectively, the remaining parameters were preserved at the default settings. The grid box was tailored to fit the active site of NRP1 protein by taking the centroid of HRG/Arg-1 as a centre of the grid, with a grid size of 15 Å × 15 Å × 15 Å and a value of spacing equal to 1 Å. The grid was centred at 32.485, –11.49, and 27.507 (in x, y, and z dimensions, respectively). The top nine scoring poses were estimated via visual investigation and the lowest docking score was selected as a representive mode.

2.4. Molecular dynamics simulations

The most effective NPASS compounds complexed with the NRP1 protein underwent molecular dynamics (MD) simulations using AMBER16 software [40]. The details of MD simulations are described elsewhere [41–45]. Briefly, the inspected NPASS compounds were defined utilizing the general AMBER force field (GAFF2) [46], while the NRP1 protein was characterized using the AMBER force field of 14SB [47].

MD simulations were conducted employing both implicit and explicit water solvents. In the implicit water solvent MD simulations, the atomic partial charges of the inspected NPASS compounds were estimated using an AM1-BCC method with the aid of the Antechamber tool within the AMBER16 [48]. Periodic boundary conditions were not utilized, and the cutoff value was set to 999 Å. Additionally, the solvent model (igb = 1) was used to present the solvation contributions [49]. All the investigated systems were first minimized for 500 steps. The minimized complexes were warmed progressively from 0 to 300 K over 10 ps. Thereafter, the NPASS-NRP1 complexes were subjected to an equilibrium stage of 50 ps, followed by a production stage of 1 ns.

In explicit water solvent MD simulations, the restricted electrostatic potential (RESP) approach was employed to determine the partial charge of each atom of the inspected NPASS compounds at the HF/6-31G* level of theory [50]. Quantum mechanics calculations were performed using Gaussian09 software [51]. TIP3P water molecules were utilized to solvate NPASS-NRP1 complexes in an octahedron box with a minimal distance of 12 Å [52]. All investigated systems were minimized by employing 5000 cycles before heating from 0 to 300 K over 50 ps under the NVT ensemble. NPASS-NRP1 complexes were then equilibrated over 1 ns under the NPT ensemble. The equilibrated coordinates of the systems were subjected to production runs of 10 and 200 ns. Atomic coordinates were recorded every 10 ps for post-dynamics analyses and binding free energy estimations. The Particle Mesh Ewald (PME) approach was applied to treat long-range electrostatic interactions [53]. Berendsen barostat with a 2 ps relaxation duration was employed to adjust the atmospheric pressure [54]. To keep temperature steady at 298 K, Langevin dynamics with gamma_ln set to 1.0 for collision frequency was applied [55]. The SHAKE algorithm was applied to constrain all bonds involving hydrogen bonds at the time step of 2 fs [56]. The GPU version of pmemd (pmemd.cuda) was applied to perform MD simulations. The molecular visualizations were generated with the help of BIOVIA Visual Studio2020 [57].

2.5. MM-GBSA binding energy

Binding energies of inspected NPASS compounds in complex with NRP1 were calculated using a molecular mechanics-generalized Born surface area (MM-GBSA) approach [58]. An updated generalized born model including Onufriev (igb = 2) was used for computing the polar solvation energy [59]. The MM-GBSA binding energies were computed based on statistically independent snapshots gathered from MD simulations in accordance with the following equation: $\mathrm{\Delta }{\mathit{G}}_{\mathbf{binding}}={\mathit{G}}_{\mathbf{complex}}-({\mathit{G}}_{\mathbf{NPASS}}+{\mathit{G}}_{\mathbf{NRP}\mathbf{1}})$

The energy term (G) was determined as: $\mathit{G}={\mathit{G}}_{\mathbf{GB}}+{\mathit{G}}_{\mathbf{SA}}+{\mathit{E}}_{\mathbf{ele}}+{\mathit{E}}_{\mathbf{vdw}}$ Through employing generalized born (GB) models, the electrostatic solvation energy (ΔG_GB) in the abovementioned equation was determined with G_SA standing for nonpolar solvation-free energy. E_ele refers to electrostatic energy. E_vdw is van der Waals energy. The entropic contributions are often disregarded because of high computational costs [60,61].

2.6. Drug-likeness characteristics

The physicochemical characteristics of the identified compound were estimated utilizing a SwissADME online website (http://www.swissadme.ch). The estimated characteristics included topological polar surface area (TPSA), relative molecular weight (MW), number of hydrogen bond donors (nOHNH), number of hydrogen bond acceptors (nON), and the partition coefficient (MlogP) [62].

2.7. ADMET characteristics

ADMET properties of the identified compounds were evaluated by an online pkCSM web server [63]. The absorption (A) was estimated by human intestinal absorption (HIA) and Caco-2 permeability. Distribution (D) was predicted based on blood–brain barrier (BBB) permeability and steady-state distribution volume (VDss). The metabolism (M) was detected using CYP3A4 inhibitor/substrate. Excretion (E) is specified according to the total clearance of the drug. Finally, the toxicity (T) was determined based on hepatotoxicity.

3. Results and discussion

3.1. In-silico protocol assessment

The performance of AutoDock Vina1.1.2 software in predicting the correct ligand-NRP1 binding pose was assessed. For assessment purposes, the co-crystallized HRG/Arg-1 inhibitor was re-docked in NRP1 binding site, and the predicted docking pose was compared to the experimental binding mode (PDB code: 5ijr [26]) (Figure 2). As depicted in Figure 2, the docking pose was similar to the native structure of the co-crystallized HRG/Arg-1 ligand with an RMSD value of 0.23 Å. AutoDock Vina results showed a docking score of –5.6 kcal/mol, forming six significant hydrogen bonds with TRP301 (2.01 Å), ASP320 (1.65, 3.11 Å), SER346 (2.60 Å), GLU348 (2.44 Å), and ILE415 (2.79 Å). Based on these results, AutoDock Vina software was utilized to screen the NPASS database for potential NRP1 inhibitors.

Figure 2. (a) (i) Superimposed structures of the co-crystalized mode (in grey) and the predicted docking pose (in green) and (ii) 2D molecular interaction of L-homoarginine (HRG/Arg-1). (b) (i) 3D and (ii) 2D representations of the predicted binding mode for kamolonol (NPC146388) within the NRP1 active site. The predicted docking score in kcal/mol is also listed.

3.2. Virtual screening of the NPASS database

Molecular docking calculations were executed via two stages of accuracy that were named fast and expensive docking stages (see computational methodology section for details). The NPASS database, containing greater than 35,000 compounds, was virtually screened using fast docking computations. The docking scores of NPASS compounds were predicted and compared to that of the co-crystallized HRG/Arg-1 ligand (calc. –5.6 kcal/mol). NPASS compounds with docking scores <−7.0 kcal/mol were selected and subjected to expensive docking computations, giving a total number of 1,332 compounds. Due to a large number of investigated compounds, a threshold value of –7.0 kcal/mol was selected to shortlist the prospective NRP1 inhibitors. Table S1 provides the evaluated fast and expensive docking scores for the top 1,332 NPASS compounds against the NRP1 protein. Figure S1 depicts 3D and 2D molecular interactions of the top four potent NPASS compounds with essential residues within the active site of the NRP1 protein. Table 1 summarizes the obtained docking scores, 2D chemical structures, and binding features of the top four potent NPASS compounds within the active site of NRP1.

Table 1. Computed fast and expensive docking scores (in kcal/mol), 2D chemical structures, and binding features for the top four potent NPASS compounds against the NRP1 protein.^a

			Docking Score (kcal/mol)
No.	NPASS Code	2D Chemical Structure	Fast	Expensive	Binding Features^b
	HRG/Arg-1^c		–5.6	–5.7	TRP301 (2.01 Å), ASP320 (1.65, 3.11 Å), SER346 (2.60 Å), GLU348 (2.44 Å), ILE415 (2.79 Å)
1	NPC146388 (Kamolonol)		–9.2	–9.3	ASN313 (2.25 Å), ASP320 (2.68 Å), LYS347 (3.01 Å), GLU348 (2.41 Å)
2	NPC108434		–8.8	–8.9	ASN313 (2.66, 2.15 Å), GLY318 (2.61 Å)
3	NPC75315 (Stelletin A)		–8.8	–8.7	SER298 (2.13 Å), ASN300 (2.26 Å), TRP301 (2.28 Å)
4	NPC19013 (Sculponin D)		–8.7	–8.7	ASP320 (2.45 Å), TYR353 (2.36 Å)

^a Data ranked based on the expensive docking scores.

^b Only hydrogen bonds are listed in Å.

^c Reference ligand.

Notably, those four NPASS compounds were filtered according to further binding energy calculations over 10 ns MD simulations. As illustrated in Table 1 and Figure S1, the four NPASS compounds and HRG/Arg-1 demonstrated approximately the same binding modes within the NRP1 active site. The most favorable interactions between the identified NPASS compounds and NRP1 active site were conventional hydrogen bonds and π-based interactions.

Kamolonol (NPC146388), isolated from F. assafoetida gum resin [64], belongs to the class of organic compounds known as coumarins. Kamolonol possesses antiplasmodial characteristics, anti-cellular hypertrophic, and anti-fibrotic effects [65,66]. Besides, kamolonol demonstrated a promising docking score towards NRP1 protein with a value of –9.3 kcal/mol (Table 1). The potency of kamolonol as an NRP1 inhibitor can be traced to its ability to exhibit multiple hydrogen bonds with the key amino acids of the NRP1 active site (Figure 2). Inspecting the binding mode of kamolonol revealed that the hydroxyl (OH) group of 2-methylcyclohexan-1-ol ring formed a hydrogen bond with the carbonyl (C = O) group of ASP320 amino acid with a bond length of 2.68 Å (Figure 2). Besides, the carbonyl (C = O) group of the methylhexanone ring exhibited three hydrogen bonds with the NH group of ASN313 and GLU348 and NH2 group of LYS347 with bond lengths of 2.25, 2.41, and 3.01 Å, respectively (Figure 2). In addition, kamolonol demonstrated π-donor hydrogen bond with TRP301 residue and π-alkyl interactions with TYR297 and TYR353 residues.

3.3. Molecular dynamics simulations

Molecular dynamics (MD) simulations of the target-inhibitor complex are essential for inspecting conformational steadiness and changes, internal motions, and the reliability of target-inhibitor binding affinities [67,68]. Consequently, MD simulations were employed herein to investigate the stability and binding affinity of the most potent NPASS compounds complexed with NRP1 protein. Considering computational cost and time, the top 5% promising NPASS compounds (i.e. 67 compounds with expensive docking scores < –8.1 kcal/mol) in complex with NRP1 were subjected to implicit water solvent MD simulations. Additionally, the corresponding MM-GBSA binding energies (ΔG_binding) were estimated and listed in Table S2. As data enrolled in Table S2, four NPASS compounds demonstrated greater binding affinity compared to the co-crystallized HRG/Arg-1 ligand (calc. –28.7 kcal/mol) with binding energies with values in the range of –21.5 to –29.8 kcal/mol. Those four potent NPASS compounds in complex with NRP1 were subjected to explicit water solvent MD simulations for 10 ns for more reliable binding affinities. The corresponding binding energies (ΔG_binding) were calculated and illustrated in Figure 3.

Figure 3. Estimated MM-GBSA binding energies for the four promising NPASS compounds against the NRP1 active site over 1 ns in implicit water solvent MD simulations and 10 and 200 ns in explicit water solvent MD simulations.

Interestingly, only kamolonol exposed binding energy lower than that of the HRG/Arg-1 (calc. –29.8 kcal/mol) over the MD course of 10 ns. A MD simulation for kamolonol was also prolonged to 200 ns, and the corresponding binding energy (ΔG_binding) was calculated (Figure 3). No discernible variation in the predicted binding energies of kamolonol within the NRP1 active site over the different simulated times of 10 and 200 ns MD simulations. Compared to HRG/Arg-1, kamolonol complexed with NRP1 demonstrated lower binding energy over 200 ns MD simulation with ΔG_binding values equal to –34.5 and –32.0 kcal/mol, respectively (Figure 3).

Inherent driving forces in the binding of kamolonol and the HRG/Arg-1 were analyzed utilizing MM-GBSA binding energy decomposition (Figure 4). Remarkably, the van der Waals forces (ΔE_vdw) demonstrated a significant contribution to the binding energy of kamolonol- and HRG/Arg-1-NRP1 complexes, with values of –38.3 and –32.5 kcal/mol, respectively (Figure 4).

Figure 4. MM-GBSA binding energy components for kamolonol (NPC146388) and the co-crystallized HRG/Arg-1 ligand with the NRP1 protein over the 200 ns MD simulations.

The electrostatic forces (ΔE_ele) were also adequate for kamolonol and HRG/Arg-1 complexed with NRP1, with average values of –21.2 and –19.6 kcal/mol, respectively (Figure 4). The abovementioned findings give statistical information on the binding energies of kamolonol as a putative NRP1 inhibitor.

To examine kamolonol-NRP1 interactions and the contribution of essential residues inside the active site of NRP1, estimated ΔG_binding values were divided into individual amino acid participations (Figure 5a). Particularly, amino acids with energy contributions less than –0.5 kcal/mol were depicted in Figure 5a. Per-residue energy decomposition analysis demonstrated that TRP301, SER346, ILE415, and ASP320 interacted with kamolonol and HRG/Arg-1 (Figure 5a). For instance, ASP320 amino acid residue participated in the total binding energy (ΔG_binding) with values of –2.6 and –2.1 kcal/mol, for kamolonol- and HRG-NRP1 complexes, respectively (Figure 5a). The average structures for kamolonol and HRG/Arg-1 inside the NRP1 active site over the 200 ns MD simulations are also displayed in Figure 5b. Kamolonol and HRG/Arg-1 formed four and two hydrogen bonds, respectively, with the key residues of the active site of NRP1 protein throughout the 200 ns MD simulations.

Figure 5. (a) The energy contribution of the most essential amino acids to the total MM-GBSA binding energy and (b) 2D representations of the anticipated binding modes for (i) kamolonol (NPC146388) and (ii) HRG/Arg-1 with the NRP1 active site in accordance with the average structure of the 200 ns MD simulations.

Inspecting the average structure for kamolonol-NRP1 complex showed that kamolonol demonstrated essential hydrogen bonds with key amino acid residues over the MD simulation. Notably, these interactions were absent in the initial predicted binding modes from docking computations. This demonstrates the importance of molecular dynamics simulations toward reliable results. Among these interactions are hydrogen bonds with ASP320, ILE415 and GLU348 (Figure 5b).

3.4. Post-dynamics analyses

3.4.1. Binding energy per-frame

To investigate the stability of kamolonol within the NRP1 active site, the correlation between the binding energy and simulation time was inspected (Figure 6). As illustrated in Figure 6, kamolonol and HRG/Arg-1 were generally stable from 32 ns until the termination of the simulation, demonstrating binding energy equal to –34.5 and –32.0 kcal/mol, respectively. These study results elucidated the steadiness of the investigated complexes throughout the 200 ns MD simulations.

Figure 6. The estimated binding energies of kamolonol (NPC146388) (in red) and HRG/Arg-1 (in black) within the NRP1 protein during the 200 ns MD simulations.

3.4.2. Hydrogen bond analysis

The hydrogen bonding of the ligands with the NRP1 protein was estimated throughout the MD course of 200 ns. The correlation between the number of hydrogen bonds and simulation time was inspected and presented in Figure S2. The most striking aspect of Figure S2 is that the average number of hydrogen bonds was 3 and 2 for the kamolonol and HRG/Arg-1. The current findings clearly supported the steady for the kamolonol and HRG/Arg-1.

3.4.3. Centre-of-mass distance

Centre-of-mass (CoM) distance was used to inspect the steadiness of the kamolonol and HRG/Arg-1 within NRP1. The CoM distance was measured between the kamolonol and HRG/Arg-1 and the ASP320 residue over the 200 ns MD simulations (Figure 7). The CoM distance was steady for kamolonol and HRG/Arg-1 complexed with NRP1 at average distances of 10.7 and 8.6 Å, respectively. The presented results demonstrated the high stability of kamolonol and HRG/Arg-1 with NRP1.

Figure 7. CoM distance (in Å) of kamolonol (NPC146388) (in red) and HRG/Arg-1 (in black), and ASP320 of NRP1 protein active site over the 200 ns MD simulations.

3.4.4. Root-mean-square deviation

To examine the structural changes of NRP1 after complexation with kamolonol and HRG/Arg-1, the root-mean-square deviation (RMSD) analysis was measured (Figure 9). The average RMSD values were 0.13 and 0.10 nm for the kamolonol- and HRG/Arg-1-NRP1 complexes, respectively, over the 200 ns MD simulations (Figure 8). The RMSD analysis showed that the two complexes preserved their stability until the end of the simulations. These results confirmed that kamolonol and HRG/Arg-1 are tightly bonded and do not disturb the overall structural constancy of NRP1 besides keeping structural integrity.

Figure 8. RMSD of the backbone atoms from the first frame of kamolonol (NPC146388) (in red) and HRG/Arg-1 (in black) with NRP1 protein during the simulation time of 200 ns.

3.5. Drug-likeness characteristics

The SwissADME website was utilized to predict the drug-likeness characteristics of kamolonol compared to HRG/Arg-1. For kamolonol the MLogP value was found to be less than five (calc. 3.0), indicating that this compound has high lipophilicity. The molecular weight of kamolonol was found to be less than 500 (calc. 398.49 g/mol). Additionally, the number of hydrogen bond acceptors (nON) was 5, and the number of hydrogen bond donors (nOHNH) was 1 (Figure 9). The TPSA of the kamolonol was 76.74, which indicates that kamolonol has good bioavailability (Figure 9). Finally, the calculated %ABS was 82.5%, explaining that kamolonol has good cell membrane permeability and oral bioavailability. Comparing the drug-likeness characteristics of kamolonol with those of HRG/Arg-1 demonstrated that kamolonol has better physiochemical properties and oral bioavailability than HRG/Arg-1.

Figure 9. Predicted physiochemical characteristics of kamolonol (NPC146388) and HRG/Arg-1 as NRP1 inhibitors.

3.6. ADMET characteristics

The pharmacokinetic properties of kamolonol were predicted and compared to HRG/Arg-1 utilizing the pkCSM online server. The pharmacokinetic properties included absorption, in which kamolonol demonstrated high Caco2 permeability (calc. 0.839) and high human intestinal absorption (HIA) (calc. 98.3), suggesting that kamolonol would be highly absorbed (Table 2). In contrast to HRG/Arg-1, HRG/Arg-1 manifested poor Caco2 permeability and HIA with values of –0.351 and 24.6, respectively (Table 2). To examine the drug distribution, the VDss value was expected with values of 0.351 and –0.438 for kamolonol and HRG/Arg-1, respectively (Table 2). The metabolism results showed that kamolonol could work as an inhibitor for the CYP3A4 enzyme; however, HRG/Arg-1 was not an inhibitor for the CYP3A4 enzyme (Table 2). The excretion was predicted based on the total clearance with values of 0.433 and 0.589 for kamolonol and HRG/Arg-1, respectively (Table 2). According to hepatotoxicity, both kamolonol and HRG/Arg-1 were not toxic (Table 2). The most striking observation to emerge from the data comparison was the superior ADMET features for kamolonol over HRG/Arg-1. These findings confirmed that the kamolonol could be used as a promising anticancer drug candidate.

Table 2. The expected ADMET properties for the kamolonol (NPC146388) and HRG/Arg-1.

	Absorption (A)		Distribution (D)	Metabolism (M)	Excretion (E)	Toxicity (T)
NPASS Code	Caco2 Permeability (cm/s)	Human Intestinal Absorption (HIA)	VDss (human)	CYP3A4 Inhibitor/Substrate	Total Clearance	Hepatotoxicity
Kamolonol (NPC146388)	0.839	98.3	0.351	Yes	0.433	No
HRG/Arg-1	–0.351	24.6	–0.438	No	0.589	No

3.7. Kamolonol vs. EG00229

EG00229, (S)−2-(3-(benzo[c][1,2,5]thiadiazole-4-sulfonamido)thiophene-2-carboxamido)−5-guanidinopentanoic acid, is an experimental inhibitor of NRP1 with IC₅₀ value of 3 μM [69]. To estimate the prospectivity of the identified NPASS compound, the binding score and mode of EG00229 towards NRP1 protein were evaluated and compared to kamolonol (NPC146388). Furthermore, MD simulations over 200 ns pursued by binding affinity estimations were performed for EG00229 complexed with NRP1 (Table 3 and Figure S3). As depicted in Figure S3, the docking pose of EG00229 within the binding pocket of NRP1 was similar to the binding mode of kamolonol, exhibiting seven hydrogen bonds with ASP320 (2.54, 2.68 Å), SER346 (2.78, 2.82 Å), THR349 (1.91 Å), and ILE415 (2.96, 3.01 Å). Compared to kamolonol (docking score = −9.3 kcal/mol), EG00229 demonstrated a good docking score against NRP1 with a value of −6.1 kcal/mol (Table 3).

Table 3. Predicted docking score (in kcal/mol) and MM/GBSA binding energies during the 200 ns MD course for kamolonol and EG00229 complexed with NRP1 protein.

		Calculated MM/GBSA binding energy (kcal/mol)
NPASS Code/Name	Docking Score (kcal/mol)	ΔEvdw	ΔEele	ΔEGB	ΔESUR	ΔGgas	ΔGSolv	ΔGbinding
Kamolonol (NPC146388)	–9.3	–38.3	–21.2	29.7	–4.8	–59.5	25.0	–34.5
EG00229	–6.1	–30.0	–18.8	23.3	–3.1	–48.9	20.2	–28.7

To gain more reliable results, EG00229 in complex with NRP1 was subjected to 200 ns MD simulations, and the corresponding MM-GBSA binding affinity was computed (Table 3). From Table 3, the average MM-GBSA binding energy (ΔG_binding) for EG00229 in complex with NRP1 was –28.7 kcal/mol, compared to –34.5 kcal/mol for the kamolonol-NRP1 complex. Similar to kamolonol, the binding affinity of EG00229 complexed with NRP1 protein was predominated by van der Waals (E_vdw) with an average value of –30.0 kcal/mol (Table 3).

4. Conclusion

In the current study, the NPASS database containing more than 35,000 natural products was virtually screened against the NRP1 protein for identifying potent inhibitors as anticancer drug candidates using the molecular docking technique. Based on docking computations, 67 compounds with expensive docking scores < –8.1 kcal/mol in complex with NRP1 were subjected to MD simulations. The computed MM-GBSA binding energy over an MD course of 200 ns demonstrated superior binding affinity of kamolonol (NPC146388) compared to the HRG/Arg-1 ligand with NRP1, demonstrating ΔG_binding values of –34.5 and –32.0 kcal/mol, respectively. During 200 ns MD simulations, energetic and structural analyses indicated perfect constancy for the kamolonol. Furthermore, kamolonol elucidated favorable physiochemical and pharmacokinetic properties compared to HRG/Arg-1. Finally, these results confirmed that the kamolonol could be used as a promising NRP1 inhibitor and hold promise for further experimental assays toward cancer treatment.

Acknowledgements

Ahmed M. Shawky would like to thank the Deanship of Scientific Research at Umm Al-Qura University for supporting this work by Grant Code: (22UQU4331174DSR43). The computational work was completed with resources supported by the Science and Technology Development Fund (STDF-Egypt, Grants No. 5480 and 7972), and Bibliotheca Alexandrina (http://hpc.bibalex.org). Dr. Mahmoud A. A. Ibrahim extends his appreciation to the Academy of Scientific Research and Technology (ASRT), Egypt, for funding Graduation Projects conducted at CompChem Lab, Egypt.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available in the Supporting Information Material of this paper.

Additional information

Funding

This work was supported by Umm Al-Qura University: [grant no 22UQU4331174DSR43].