Computational study of novel natural inhibitors targeting aminopeptidase N(CD13)

Virtual screening of natural products database against CD13

The S1 pocket is an important regulatory site of CD13, as the Val (P1) side chain in the Ang IV complex binding to this active site could cause vasodilatation, hypertrophy, and stimulate NF-κB pathway, and then neovascularization occurs. Therefore, this pocket region was selected as a reference site. A total of 17,931 biogenic-for sale-named product molecules were taken from the ZINC15 database. Chemical structure of the receptor protein, CD13, was selected to test the pharmacologic properties of other compounds. One inhibitor, Bestatin, was chosen as the reference compound. After screening, 4761 compounds were found to be eligible to bind stably with CD13; among those, 318 compounds were found to have higher LibDock scores than Bestatin. (LibDock score: 97.2164), which are listed in Supplementary Table 1. The top 20 ranked compounds are listed in Table 1.

ADME and toxicity prediction

Pharmacologic properties of all the 20 selected ligands with Bestatin were predicted using the ADME module of Discovery Studio 4.5, including aqueous solubility level, blood-brain barrier level, CYP2D6 binding, hepatotoxicity, human intestinal absorption level, and plasma protein binding properties (Table 2). The aqueous solubility prediction (defined in water at 25°C) indicated that all the compounds were soluble in water. Only one compound, ZINC000018185774, was predicted to be inhibitor with CYP2D6, which was an essential enzyme in drug metabolism. As to hepatotoxicity, all compounds were found to be toxic, while Bestatin is nontoxic. For human intestinal absorption, except 2 compounds, ZINC000006092198 and ZINC000018169010, other ligands had as absorption level as good as Bestatin did. Plasma protein binding properties indicated that ZINC000006094124, ZINC000014820552, ZINC000006092198, ZINC000000899884, ZINC000014820405, ZINC000020470300, ZINC000006092239, ZINC000014820583 and ZINC00000 5732375 had good absorption, while other 11 compounds had weak absorption.

Safety and toxicity should also need full consideration in this study. To confirm the safety of the top20 compounds, different kinds of toxicity indicators of these compounds and Bestatin, including rodent carcinogenicity (based on the U.S. National Toxicology Program dataset), developmental toxicity potential properties, Ames mutagenicity were predicted by using a computational method in the TOPKAT module of Discovery Studio 4.5 (Table 3). Results indicated that 13 compounds were found to be nonmutagenic, and 2 compounds showed not to have developmental toxicity potential(DTP), while one compound, ZINC000000895551, had little DTP. The reference Bestatin was predicted to have low rodent carcinogenicity and low DTP. Considering all the results in Table 3, ZINC000000895551 and ZINC000014820583 were identified to be the ideal leading compounds with non-CYP2D6 inhibitors and less rodent carcinogenicity, together with Ames mutagenicity, and developmental toxicity potential compared with other compounds. Therefore, ZINC000000895551 and ZINC000014820583 were confirmed as safe drug candidates and were selected for the subsequent research (Figure 1).

Figure 1. Structures of novel compounds selected from virtual screening.

Analysis of ligand binding and ligand pharmacophore

To identify ligand binding mechanisms of these compounds and Bestatin, ZINC000000895551 and ZINC000014820583 were docked into the molecule structure of CD13 using CDOCKER module, and CDOCKER potential energy was calculated and displayed in Table 4. RMSD between the docked pose and the crystal structure of the Bestatin and CD13 complex was 0.75 Å, proving that the CDOCKER module applied in this study was strongly reliable for reproducing the experiment. As results shown in Table 5, the CDOCKER potential energy of ZINC000000895551 and ZINC000014820583 were significantly lower than the reference ligand Bestatin (ZINC000001532730) (24.7496 kcal/mol), which suggested that these 2 compounds had a much higher binding affinity with CD13 than Bestatin. Hydrogen bonds and π-related interactions were also performed using structural computation study (Figures 2 and 3). Results showed that ZINC000000895551 formed 6 pairs of hydrogen bonds with CD13, by the H²⁸ of the compound with ALA819:O of CD13, H³² of the compound with GLU364:OE2 of CD13, the O²⁰ of the compound with ARG305:HD1, the H²⁵ of the compound with GLU364:OE1, the H³⁰ of the compound with SER307:OG and the H³¹ of the compound with THR820:O. Besides, 4 π-related interactions were formed in the complex, including the π-Orbitals between ZINC000000895551 and the ALA308, LEU368, LEU821, and VAL822 of CD13. ZINC000014820583 formed 6 pairs of hydrogen bonds with CD13, by the H³⁴ of the compound with ALA819:O of CD13, the O⁹ of the compound with ASN823:HN of CD13, the H²⁴ of the compound with LEU368:O of CD13, the H³⁰ of the compound with GLU364:O of CD13, the H³⁶ of the compound with SER307:OG of CD13 and the H³⁷ of the compound with THR820:O of CD13. It also formed 5 pairs π-related interactions in the complex, consist of 2 pairs between ALA308 of CD13 and the compound, ARG305 of CD13 with the compound, LEU368 of CD13 with the compound, and LEU821 of CD13 with the compound. As to the reference compound Bestatin, it formed 2 hydrogen bonds with CD13, by the H³⁸ of the compound with THR820:OG1 of CD13, and H³⁶ of the compound with THR820:O of CD13. Also three π-related interactions formed between Bestatin and CD13, including ALA308, LEU368 and LEU821 of CD13 with Bestatin.

Figure 2. Schematic of intermolecular interaction of the predicted binding modes of (A) ZINC000014820583 with aminopeptidase N(CD13), (B) ZINC000000895551 with CD13.

Figure 3. Schematic drawing of interactions between ligands and aminopeptidase N(CD13). Blue represents positive charge; red represents negative charge; and ligands are shown in sticks, with the structure around the ligand-receptor junction shown in thinner sticks. (A) ZINC000014820583-CD13 complex. (B) ZINC000000895551-CD13 complex.

Molecular dynamics simulation

In order to estimate the stabilities of the ligand-CD13 complexes on natural environment conditions, the molecular dynamics simulation module test was conducted. According to the original conformations from the molecular docking experiment through CDOCKER module, RMSD curves and potential energy chart of each complex were shown in Figure 4. The trajectories of all complexes achieved equilibrium after 18 ps; RMSD and potential energy of these complexes tends to stability with time. The molecular dynamics simulation result suggested that these hydrogen bonds and π-related interactions formed by compounds and CD13 contributed to the stability of these complexes. Results showed that ZINC000000895551 and ZINC000014820583 could interact with CD13, and their complexes with CD13 could exist in a natural environment steadily, and act as regulatory roles to CD13 just like the reference Bestatin did.

Figure 4. Results of molecular dynamics simulation of the compounds ZINC000014820583 and ZINC000000895551. (A) Potential energy. RMSD, root-mean-square deviation. (B) Average backbone root-mean-square deviation.

Discovery studio software and ligand library

Discovery Studio 4.5 software (BIOVIA, San Diego, California, USA) is a suite of software for simulating small molecule and macromolecule systems; it is a new generation of molecular modeling and environmental simulation software for the life sciences field, which aims to provide protein modeling, optimization, and drug design tools by applying protein structure and structural biologic computation. Numerous leading compounds and drug candidates were identified and refined through this method. LibDock module of Discovery Studio was employed for virtual screening; CDOCKER module was used for docking study; and ADME module was analyzed for pharmacologic properties. The Natural Products database in the ZINC15 database was selected to screen STING agonists. The ZINC15 database is a free database of commercially available compounds provided by the Irwin and Shoichet Laboratories, Department of Pharmaceutical Chemistry, University of California, San Francisco (San Francisco, California, USA).

Structure-based virtual screening using LibDock

Ligand-binding pocket region of CD13 was selected as the binding site to screen potential compounds to inhibit CD13. Virtual screening was carried out using the LibDock module of Discovery Studio 4.5 [21]. LibDock is a rigid-based docking module. It calculates hotspots for the protein through a grid placed into the binding site and polar and apolar probes. The hotspots are further used to align the ligands to form favorable interactions. The Smart Minimiser algorithm and CHARMM force field (Harvard University, Cambridge, Massachusetts, USA) were performed for ligand minimization, and all ligand poses were ranked based on the ligands score. The 1.9-Å crystal structure of human CD13 (Protein Data Bank identifier: 4fyq) was downloaded from the Protein Data Bank and the inhibitor Bestatin (Zinc15 database identifier: zinc000001532730) was downloaded from Zinc15 database, and they were imported to the working circumstance of LibDock. The chemical structure of CD13 is shown in Figure 5. The protein was prepared by serval steps including removing crystal water and other heteroatoms around it, followed by addition of hydrogen, protonation, ionization, and energy minimization. The CHARMM force field and the Smart Minimiser algorithm were applied for energy minimization [22]. The minimization performed 2000 steps with a root mean square gradient tolerance of 12.277, and the final root mean square gradient was 0.690. The prepared protein was prepared to define the binding site. Using the ligands (Bestatin) binding position, the active site for docking was generated. Virtual screening was performed by docking all the prepared ligands at the defined active site using LibDock. Based on the LibDock score, all the docked poses were ranked and grouped, and all compounds were ranked according to the LibDock score.

Figure 5. Molecular structure of aminopeptidase N(CD13).

Absorption, distribution, metabolism, and excretion and toxicity prediction

The ADME module of Discovery Studio 4.5 was employed to calculate absorption, distribution, metabolism, and excretion (ADME) of selected compounds, including their aqueous solubility, blood-brain barrier penetration, cytochrome P-450 2D6 (CYP2D6) inhibition, hepatotoxicity, human intestinal absorption and plasma protein binding level. TOPKAT module of Discovery Studio 4.5 was employed to calculate the toxicity and other properties of all the potential compounds, such as U.S. National Toxicology Program rodent carcinogenicity, Ames mutagenicity, developmental toxicity potential, and rat oral median lethal dose (LD50) and chronic oral lowest observed adverse effect level (LOAEL). These pharmacologic properties were fully considered in screening proper probable drug candidates for CD13.

Molecule docking and pharmacophore prediction

CDOCKER module of Discovery Studio 4.5 was used for molecular docking study. CDOCKER is a molecular docking method based on CHARMM force field, which can produce high-precision docking results. The CHARMM force field was used for both receptors and ligands. Receptor is held rigid, while ligands are allowed to flex during the docking process. For each complex pose, the CHARMM energy and the interaction energy, which indicated ligand binding affinity, were calculated. Crystal structure of CD13 was obtained from the protein data bank. The crystal water molecules were generally removed in a rigid and semi-flexible docking process, causing the fixed water molecules to possibly affect the conformation of the receptor-ligand complex [23, 24]. The water molecules were removed, and hydrogen atoms were added to the protein. To prove the reliability of the analysis result, the initial inhibitor compound Bestatin was extracted from the Zinc15 database, the same to the Natural Products screened, and then Bestatin was docked into the crystal structure of CD13 to compare the root-mean-square deviation (RMSD) with these 2 conformations. The binding site sphere of CD13 was defined as the regions coming within 5-Å radius from the geometric centroid of the ligand Bestatin. During the docking process, the ligands were allowed to bind with the residues of protein groups within the binding site sphere. The structures of identified hits were prepared and docked into the binding site of CD13. Different poses of each ligand-CD13 complex were generated and analyzed on the basis of CDOCKER interaction energy.

Molecular dynamics simulation

The best binding conformations of the ligand-CD13 complexes among the poses predicted by the molecule docking program were selected and prepared for molecular dynamics simulation. The ligand-receptor complex was put into an orthorhombic box and solvated with an explicit periodic boundary solvation water model. To simulate the physiologic environment, sodium chloride was added to the system with the ionic strength of 0.145. Then the system was subjected to the CHARMM force field and relaxed by energy minimization (500 steps of steepest descent and 500 steps of conjugated gradient), with the final root mean square gradient of 0.227. The system was slowly driven from an initial temperature of 296 K to the target temperature of 302 K for 2 ps, and equilibration simulations were run for 5 ps. Molecular dynamics simulation (production module) was performed for 25 ps with time step of 1 fs. The simulation was performed with the normal pressure and temperature system at a constant temperature of nearly 300 K during the process. The particle mesh Ewald algorithm was used to calculate long-range electrostatics, and the linear constraint solver algorithm was adapted to fix all bonds involving hydrogen. With initial complex setting as a reference, a trajectory was determined for RMSD, potential energy, and structural characteristics through the Discovery Studio 4.5 analysis trajectory protocol.