- 1Covenant University Bioinformatics Research (CUBRe), Covenant University, Ota, Ogun State, Nigeria
- 2Department of Chemistry, Covenant University, Ota, Ogun State, Nigeria
- 3Department of Biochemistry, Covenant University, Ota, Ogun State, Nigeria
- 4Center for Global Infectious Disease Research, Seattle Children’s Research Institute, Seattle, WA, United States
- 5Department of Computer and Information Science, Covenant University, Ota, Ogun State, Nigeria
- 6Covenant Applied Informatics and Communications ACE (CApIC-ACE), Covenant University, Ota, Ogun State, Nigeria
- 7Division of Applied Bioinformatics, German Cancer Research Center (DKFZ), Heidelberg, Germany
Plasmodium falciparum (Pf) 5-aminolevulinic acid synthase (5-ALAS) is an essential enzyme with high selectivity during liver stage development, signifying its potential as a prophylactic antimalarial drug target. The aim of this study was to identify important potential lead compounds which can serve as inhibitors of Pf 5-ALAS using pharmacophore modeling, virtual screening, qualitative structural assessment, in silico ADMET (Absorption, Distribution, Metabolism, Excretion and Toxicity) evaluation and molecular dynamics simulation. The best model of the tertiary structure of Pf 5-ALAS was obtained using MolProbity, while the following databases were explored for the pharmacophore-based virtual screening: CHEMBL, ChemDiv, ChemSpace, MCULE, MCULE-ULTIMATE, MolPort, NCI Open Chemical Repository, LabNetwork and ZINC databases. 2,621 compounds were screened against the modeled Pf 5-ALAS using AutoDock vina. The post-screening analysis was carried out using Discovery Studio while molecular dynamics simulation was performed on the best hits using NAMD-VMD and Galaxy Europe platform. Compound CSMS00081585868 was observed as the best hit with a binding affinity of -9.9 kcal/mol and predicted Ki of 52.10 nM, engaging in seven hydrogen bonds with the target’s active site amino acid residues. The in silico ADMET prediction showed that all ten best hits possessed relatively good pharmacokinetic properties. The qualitative structural assessment of the best hit, CSMS00081585868, revealed that the presence of two pyridine scaffolds bearing hydroxy and fluorine groups linked by a pyrrolidine scaffold contributed significantly to its ability to have a strong binding affinity with the receptor. The best hit also showed stability in the active site of Pf 5-ALAS as confirmed from the RMSD obtained during the MD simulation.
1. Introduction
Malaria, a disease caused by eukaryotic Plasmodium parasites, continues to kill upward of 6,27,000 people per year and with a further 241 million clinical cases (1). Parasite resistance to antimalarial drugs continues to pose a threat to malaria prevention and control, placing an increasing responsibility on the development of alternative therapies (2). Amongst the metabolic processes critical to Plasmodium falciparum (Pf) growth during the liver stage of development is the heme biosynthetic pathway (3–5). The rate-limiting enzyme in heme biosynthesis, 5-aminolevulinic acid synthase (5-ALAS), is responsible for catalyzing a condensation reaction between succinyl-CoA and glycine to yield 5-aminolevulinic acid (ALA) (6). Using computational network studies, Bazzani et al. (3) identified Pf 5-ALAS as an essential enzyme with high selectivity during liver stage development, signifying its potential as a prophylactic antimalarial drug target. The essentiality of the pathway is likely due to the incorporation of heme into cytochromes, part of the Pf mitochondrial electron transport chain (7).
Even though heme-dependent electron transfer is vital for Pf asexual blood-stage replication, parasite heme biosynthesis appears not to be essential. This is likely due to parasites sequestering heme from host erythrocytes to fulfill the obligatory demand of heme for viable cytochromes b, c and c1 of the electron transport chain (8). This finding and reports that parasites lacking the 5-ALAS gene can complete asexual blood stage replication show that de novo biosynthesis of heme is not essential during this life cycle stage (5). Inhibiting 5-ALAS, however, has been shown to strongly inhibit the progression of liver stage-to-blood stage transition and also prevent mosquito stage sporozoite maturation (9), and thus 5-ALAS has been suggested as a potential target for malaria prophylaxis and preventing malaria transmission (8).
Structure-based virtual screening (SBVS) presents a fast and cost-effective alternative approach to high throughput screening (HTS) in the drug development process (10). This technique can also improve the efficiency of lead compound selection by giving priority to compounds with a higher probability of success during the experimental testing phase (11). By combining pharmacophore modeling with molecular docking and in silico toxicity testing, computational approaches have proven to be potent in identifying potential drug compounds that may serve as leads in drug development (12–14). In this study, we aimed to identify important potential lead compounds which can serve as inhibitors of Pf 5-ALAS using pharmacophore modeling, virtual screening, qualitative structural assessment, in silico ADMET (Absorption, Distribution, Metabolism, Excretion and Toxicity) evaluation and Molecular Dynamics (MD) simulation.
2. Materials and methods
2.1. Plasmodium falciparum 5-ALAS structure prediction
The homology modeled 3D structure of Pf 5-ALAS was built via SWISS-MODEL using the structure of ALAS from Saccharomyces cerevisiae non-covalently bound to PLP cofactor (PDB ID: 5TXR), due to the absence of its experimental structure at the time of the research (15, 16). Ab initio models were also retrieved from AlphaFold (ID: AF-Q8I4 × 1-F1) and Robetta (17). Robetta is an automated online platform for both protein structure and function predictions based on the amino acid sequence of the protein. The amino acid sequence was retrieved from the Protein Database on the National Center for Biotechnology Information (NCBI) based with accession ID: XP_001350846.
2.2. Assessment of modeled 5-ALAS structures
Preliminary structure assessments of the predicted models of the protein were carried out using SWISS-MODEL (MolProbity and Clash scores) (18, 19), then ERRAT and VERIFY from the UCLA-DOE LAB – SAVES v6. The lower values of the MolProbity and clash score with higher percentiles in both cases indicate a better model. For the ERRAT score, an average overall quality factor of approximately 91% is a good model, while for VERIFY, at least 80% of the amino acids should score = 0.2 in the 3D/1D profile. The best model from the five models generated from Robetta was identified, while a further assessment was carried out on the best model and the AlphaFold model using MolProbity, http://molprobity.biochem.duke.edu/index.php. MolProbity is a web server employed in the validation of the quality of modeled protein and nucleic acid structures. The MolProbity score is a combined single indicator of model quality based on the clash, rotamer and Ramachandran parameters.
2.3. Alignment of the modeled Pf 5-ALAS structures
The alignment of the modeled 3D structures of Pf 5-ALAS from SWISS-MODEL, Robetta and AlphaFold was carried out using PyMOL. This was performed to measure the root mean square deviation (RMSD) between the three models using the best model from Robetta (due to the structure assessment) as the reference protein. RMSD indicates how close or similar the protein structure are. Normally, more similar proteins usually have small values of RMSD. RMSD is usually reported in Angstrom (Å).
2.4. Active site prediction of the 5-ALAS modeled structure
The active sites of the modeled 5-ALAS structure with the best structure assessment were predicted using Computed Atlas of Surface Topography of proteins (CASTp) 3.0 (20, 21) and the Prankweb server (22, 23). The modeled protein structure was submitted on the webservers. The necessary amino acids for binding interactions predicted by the two servers were compared to determine the similarity between the two predicted active sites.
2.5. Pharmacophore-based screening and ligand-library generation
A pharmacophore is a chemical framework that carries the structural features of a known active compound. A pharmacophore model is developed based on the structural features of an active compound and then employed in filtering, evaluating and screening databases of molecules (24). A ligand-based pharmacophore model based on the modeled structure of Pf 5-ALAS was developed using the Pharmit server1 (25, 26). The Pharmit server offers a set-up for the virtual screening of databases utilizing suitable pharmacophore features. The modeled Pf 5-ALAS was inputted as the receptor while the pharmacophore features were set based on the interaction of pyridoxal 5′-phosphate, a reported cofactor of Pf 5-ALAS (PubChem ID: 1051) in the binding site of the protein.
Four key properties of ligands in terms of hydrogen bond acceptors, hydrogen bond donors, hydrophobicity, and aromaticity were chosen to construct an effective pharmacophore query for virtual screening. The hit screening parameters were also set based on Lipinski’s rule of five (27): molecular weight ≤ 500, hydrogen bond acceptors ≤ 10, hydrogen bond donors ≤ 5, a logP (octanol-water partition coefficient) value ≤ 5 and Veber’s filter: rotatable bonds ≤ 10, Aromatics ≤ 2 and total polar surface area (TPSA) ≤ 140 Å. The following databases were explored for the pharmacophore-based virtual screening: CHEMBL, ChemDiv, ChemSpace, MCULE, MCULE-ULTIMATE, MolPort, NCI Open Chemical Repository, LabNetwork and ZINC databases.
2.6. Ligand preparation
The 2,755 compounds obtained from the pharmacophore-based screening were downloaded in their.sdf format and converted to their corresponding 3D structures using the OpenBabel (28) panel of PyRx (29) software. 2,621 compounds were successfully converted to the AutoDock docking format (.pdbqt) due to the inability to set up a force field for the remaining 134 compounds. Pyridoxal 5′-phosphate was also added to the ligand library of 2,621 compounds for the docking simulation.
2.7. Active site determination
The 5-ALAS active site was predicted using CastP and PrankWeb online servers. The active site residues generated from the prediction were used to construct a grid box around the active site using 100, 100, and 100 as the number of points in x, y, and z directions. A grid point spacing of 0.375 Å was employed, with the center grid box set at 31.369, 26.246, and 6.703 Å of the modeled Pf 5-ALAS structure.
2.8. Virtual screening and post-screening analyses
The molecular docking analysis was performed on Autodock vina installed on the Covenant University Bioinformatics Research (CUBRe) high-performance computer (HPC). The post docking analysis was carried out using the Discovery Studio 2021 Client (30). The structural representations of the best hits in relation to the binding affinities were also examined to build the SAR (Structure-Activity Relationship) studies. SAR refers to approaches built to find relationships between chemical structures of studied compounds and biological activity (or target property). The extensive qualitative structural assessment was carried out on the four best hits to examine the heterocyclic templates and functional groups responsible for the different binding interactions observed in the post-screening analysis.
2.9. In silico ADMET prediction
The increased drug attrition rate has been attributed to toxicity and poor pharmacokinetics (31). ADMET prediction is employed in estimating the pharmacokinetic properties and toxicity risk (32). This also predicts if proposed lead compounds stand the chance of being orally active drugs. The OSIRIS Property Explorer tool was employed in this study for the prediction of ADMET of the best hits. The following pharmacokinetic properties were examined: molecular weight, solubility (log S), hydrophilicity (log P), topological polar surface area (TPSA), drug-likeness, and drug score. Furthermore, the toxicity risks of the compounds in the form of tumorigenic, mutagenic, irritating, and reproductive risks were also examined.
2.10. Molecular dynamics simulation
The MD simulation was carried out by using NAMD (Nanoscale Molecular Dynamics) (33) and VMD (Visual Molecular Dynamics) (34). This helped to examine the stability of the best confirmation of CSMS00081585868 from the docking studies in the binding site of the modeled Pf 5-ALAS. Energy minimization of 1,000 steps and production run of 5,00,000 steps (1 ns) were employed in the simulation. The topology of the target was generated using VMD while that of ligand built using Charmm36 forcefield of the Charmm-GUI webserver. The topologies (pdb and psf) of the Pf 5-ALAS and CSMS00081585868 were merged, and the complex was solvated using VMD. The simulation was performed at constant pressure of 1 atm and temperature of 310 K using Periodic Boundary conditions. All the necessary parameters for the simulation were defined in a script and executed using NAMD. The RMSD (root mean square deviation), RMSF (root mean square fluctuation) and PCA (principal component analysis) of the simulation results were carried out using Bio3D on the Galaxy Europe platform (35) while the hydrogen bond (h-bond) analysis was carried out using VMD software.
3. Results
3.1. Assessment of modeled 5-ALAS structure
Structure assessment of the Pf 5-ALAS models obtained from SWISS-MODEL, Robetta and AlphaFold revealed that the Robetta modeled structure was the best based on the different assessment scores. Robetta generated five models of the protein with a confidence score of 0.67. Modeled structure (Robetta_3) was selected as the best model based on the parameters obtained from the structure assessments. Robetta_3 had a MolProbity score of 1.47 (96th percentile), clash score of 3.93 (96th percentile), ERRAT score of 94.53%, and VERIFY score of 82.70% (Table 1). However, further structure assessment of the modeled Robetta_3 in comparison with the AlphaFold model showed that modeled Robetta_3 is better (Table 2). Robetta_3 had 95.86% of the amino acids in the Ramachandran favored region and 0.64% as Ramachandran outliers, while the AlphaFold model had 93.23% of the amino acids in the Ramachandran favored region and 5.25% as Ramachandran outliers.
Table 1. Comparison of the structure assessments of the modeled Pf 5-ALAS 3D structures based on SWISS-MODEL, ERRAT and VERIFY.
Table 2. Further structure assessment evaluation of the AlphaFold generated structure and the best model from Robetta webserver using MolProbity.
3.2. Alignment of the modeled Pf 5-ALAS structures
The RMSD of the AlphaFold modeled structure was obtained to be 12.894, having 3,978 atoms different from that of reference model. Furthermore, SWISS-MODEL structure had RMSD value of 2.053 and 3,200 atoms different from the reference model. This suggested that the homology modeled structure from SWISS-MODEL was more similar to that of the reference Robetta model than that of the AlphaFold model (Figure 1).
Figure 1. Alignment of the modeled Pf 5-ALAS Structures from Robetta (green colored), SWISS-MODEL (purple colored) and AlphaFold (blue colored).
3.3. Pharmacophore-based screening and ligand-library generation
The pharmacophore features of pyridoxal 5′-phosphate gave rise to the generation of the parameters used in the pharmacophore-based modeling on the Pharmit server. The x, y, and z coordinates of the different features (Hydrogen donor, Hydrogen acceptor, Aromatics and Hydrophobic) are shown in Table 3. The pharmacophore-based screening of the different databases (CHEMBL30, ChemDiv, ChemSpace, MCULE, MCULE-ULTIMATE, MolPort, NCI Open Chemical Repository, LabNetwork and ZINC) gave rise to 2,755 hits (Table 4). The total number of hits obtained from the databases are in the order: ChemSpace > ZINC > MCULE > MolPort > CHEMBL30 > ChemDiv > LabNetwork > NCI Open Chemical Repository > MCULE-ULTIMATE.
3.4. Virtual screening analyses
A total of 2,621 compounds, including pyridoxal 5′-phosphate, were docked into the predicted active site of Pf 5-ALAS, and the best hits were obtained (Table 5). The binding affinities of the top ten hits from the virtual screening were between -9.9 and -9.1 kcal/mol, while the predicted inhibition constant (Ki) was within the ranges of 52.10 and 202.01 nM. All the top ten hits had better binding affinities than the cofactor pyridoxal 5′-phosphate, with a docking score of -6.4 kcal/mol. However, compound CSMS00081585868, obtained from the Chemspace database, had the best binding affinity with a docking score of -9.9 kcal/mol and Ki of 52.10 nM. This means that this compound has the highest probability of causing inhibition of Pf 5-ALAS, thereby limiting the activity of the parasite at the liver stage.
Table 5. The structures, binding affinities and predicted Ki of the ten best hits from the virtual screening alongside the naturally occurring cofactor pyridoxal 5′-phosphate.
3.5. Post-screening analyses
The interactions of the atoms of the best hits with the amino acid residues in the active site of Pf 5-ALAS were analyzed using Discovery Studio 2021 client. The interactions observed from the post-screening analyses include conventional hydrogen bond, carbon-hydrogen bond, pi-cation, pi-anion, pi-sigma, pi-pi stacked, alkyl and pi-alkyl. These interactions and bond lengths are presented in Table 6. All these interactions contribute significantly to the binding affinities and docking score obtained for each of the best hits. The presence of intermolecular hydrogen bonds formed between the ligand and the amino acids in the active sites contributes significantly to the strength of the complex formed, which invariably impacts the docking result positively (36). Hence, the availability of hydrogen bond acceptors (HBA) and hydrogen bond donors (HBD) in the structure of the ligand is important. According to Lipinski’s rule of five (37), a drug candidate should possess HBA ≤ 10 and HBD ≤ 5. CSMS00081585868, the best hit from the screening, has eight HBA and three HBD, which were responsible for seven hydrogen bonds between the atoms of the compound and the amino acid residues in the Pf 5-ALAS active site. The hydrogen bonds were formed at SER 43, ASP 61, LYS 63, ASN 331, TYR 328, ASP 430, and HIS 433, and these interactions contributed to the high binding affinity observed for the compound. Likewise, the atoms of MolPort-047-716-699 formed different interactions with the amino acid residues in the Pf 5-ALAS active site, which include carbon-hydrogen bond at SER 234, amide-pi stacked/pi-pi stacked interaction at ASN 194 and TYR 202, alkyl/pi-alkyl interaction at PRO 212, LYS196 and TYR 202. Similarly, CSCS00058497692, with five HBA and two HBD, formed three hydrogen bonds with LEU 462, GLY 470 and SER 463 in the Pf 5-ALAS active site. More interactions formed include pi-donor hydrogen bond at SER 326, pi-cation/pi-anion interaction with LYS 63 and ASP 44. Furthermore, pi-pi stacked interaction was also observed at TYR 328, while alkyl/pi-Alky interaction was observed at LEU 469. Furthermore, MolPort-047-733-181 formed three hydrogen bonds with the Pf 5-ALAS active site at SER 40, SER 43 and ASN 62. Pi-cation/pi-anion interactions were also observed at LYS 65 and ASP 44 while pi-pi stacked interaction was observed at TYR 328 and alkyl/pi-alkyl interactions at VAL 48, ILE 68, PRO 47, TYR 328.
Table 6. The interactions and bond lengths of the ten best hits in the active site alongside the naturally occurring cofactor pyridoxal 5′-phosphate.
3.6. Qualitative structural assessment of the four best hits in the docking model
The structural representations of the four best hits from the virtual screening were further examined to understand their structure-activity relationship (SAR). This was carried out to understand the significant heterocyclic templates and functional groups present on the compounds that contributed to the interactions formed between the atoms of the compounds and the active site of the Pf 5-ALAS. With the understanding of the SAR of the best hits, new scaffolds can be designed that can stand a chance as new potential antimalarial prophylaxis drugs against Pf 5-ALAS. Also, functional group interconversion can be employed with careful hit-to-lead optimization. The qualitative structural assessments of the best five hits from the virtual screening are discussed below:
3.6.1. CSMS00081585868
The best hit, CSMS00081585868, comprises two pyrazine rings linked together by a pyrrolidine template (Figure 2A). Pyrrolidine (also known as tetrahydropyrrole) is a nitrogen-containing five-membered heterocyclic compound, while pyridine is a six-membered ring nitrogen-containing heterocyclic compound. The two pyridine rings possess hydroxyl (OH) on position-2, an amide linkage on position-3 and fluorine (F) on position-5. The O of the OH on one of the pyrazine rings formed a hydrogen bond with ASN 331 at a bond length of 2.91 Å, while the H of the NH on both pyrazine rings formed hydrogen bonds with ASP 430 and SER 43 at bond lengths of 2.20 and 2.39 Å, respectively. The presence of the fluorine atoms on the pyrazine rings also contributed significantly to the binding affinities forming strong hydrogen bonds with SER 43 and HIS 433 at bond lengths of 2.39 and 2.59 Å. The two O of the amide functional groups on the compound formed hydrogen bonds with LYS 63 and TYR 328 at bond lengths of 5.86 and 3.09 Å (Figure 2B).
Figure 2. (A) Qualitative structural assessment of CSMS00081585868. (B) Post docking two-dimensional visualization of CSMS00081585868 from the virtual screening in the binding pockets of Pf 5-ALAS.
3.6.2. MolPort-047-716-699
The major heterocyclic templates in the structure of MolPort-047-716-699 are pyrimido-azepine and pyrrolidine, linked by a ketone (C=O) functional group (Figure 3A). The pyrimido-azepine is a fusion of pyrimidine (a six-membered ring heterocyclic compound with nitrogen on position-1 and 3) and hexahydro-1H-azepine (a saturated seven-membered ring nitrogen-containing heterocyclic compound. The O of the C=O linkage formed a carbon-hydrogen bond with SER 234 at a bond length of 3.08 Å. A van der waals interaction was also observed between the compound and LYS 195 (Figure 3B).
Figure 3. (A) Qualitative structural assessment of MolPort-047-716-699. (B) Post docking two-dimensional visualization of MolPort-047-716-699 from the virtual screening in the binding pockets of Pf 5-ALAS.
3.6.3. CSCS00058497692
CSCS00058497692, comprises two heterocyclic rings: dihydroquinoline (fusion of benzene and piperidine rings) and 1,2,3-triazole (five-membered ring with nitrogen at position-1, 2 and 3). These rings are linked together by a cyclohexane template (Figure 4A). The H of the CH2OH on position-4 of the triazole is responsible for the hydrogen bond formed with SER 463 at a bond length of 2.22 Å. The H of the N on position-3 of the triazole is responsible for the hydrogen bonds formed with LEU 462 and GLY 470 at 2.93 and 1.86 Å, respectively (Figure 4B).
Figure 4. (A) Qualitative structural assessment of CSCS00058497692. (B) Post docking two-dimensional visualization of CSCS00058497692 from the virtual screening in the binding pockets of Pf 5-ALAS.
3.6.4. MolPort-047-733-181
MolPort-047-733-181, comprises two heterocyclic rings: quinoline (fusion of benzene and pyridine rings) and oxadiazabicyclo[3.3.2]decane. These rings are linked together by a carbonyl group (C=O; Figure 5A). The H of the NH of quinoline formed a hydrogen bond with SER 43 at a bond length of 2.16 Å, while the O of the OH on position-4 of the quinolone formed a hydrogen bond with ASN 62A at 2.31 Å. The O of the oxadiazabicyclo[3.3.2]decane is responsible for the carbon-hydrogen bond formed with SER 40 at a bond length of 3.70 Å. Furthermore, the methyl group on position-7 of the quinoline template formed alkyl/pi-alkyl interactions with VAL 48, ILE 68, and TYR 328 at bond lengths of 4.99, 4.47, and 5.13 Å, respectively (Figure 5B).
Figure 5. (A) Qualitative structural assessment of MolPort-047-733-181. (B) Post docking two-dimensional visualization of MolPort-047-733-181 from the virtual screening in the binding pockets of Pf 5-ALAS.
3.7. In silico ADMET prediction
The pharmacokinetic properties and toxicity risks of all the best hits, as predicted by OSIRIS Property Explorer, are presented in Table 7. Molecular weights (MW) of compounds determine their ease of distribution within cells. Compounds with lesser MW tend to be more easily distributed than compounds with higher MW, hence the benchmark of 500 g/mol. All the best hits had MWs within the acceptable range (between 341.0 and 437.5 g/mol). The logarithm of the partition coefficient between n-octanol and water is used to determine the clog P, which is a measure of a compound’s hydrophilicity or lipophilicity. A clog P of more than 5.0 indicates low hydrophilicity or poor absorption, but clogP values less than 5.0 are acceptable. All the best hits possessed clog P values less than 5.0, which suggests good absorption potentials for all the compounds. The topological polar surface area (TPSA) of a molecule is strongly related to its hydrogen bonding and is a good predictor of bioavailability (38). A score of less than 160 Å2 is considered acceptable for TPSA, indicating that the molecule will have good oral bioavailability (39). All the best hits had acceptable TPSA scores.
Table 7. Pharmacokinetic properties and toxicity risks of the 10 best hits as estimated using the OSIRIS property explorer tool.
Solubility is also a significant consideration in pharmacokinetics as it impacts both absorption and distribution. It is measured as a logarithm of the solubility measured in mol/dm3. All the compounds possess estimated log S values greater than -4 except for ZINC5276997, and this corresponds to the score of more than 80% of marketed drugs. The drug-likeness of the compounds is predicted on a positive or negative basis. A positive number indicates that the molecule has a high percentage of fragments that are often seen in marketed drugs. All the ten best hits had positive values except for CSCS00058497692 and MolPort-045-962-768 with -3.59 and -3.32, respectively. The drug score parameter combines drug-likeness, cLogP, logS, molecular weight, and toxicity risk into one easy-to-understand number that may be used to assess a compound’s overall potential to become a drug. The higher the drug score value, the higher the compound’s chance of being considered a drug candidate (40). Compound MolPort-047-716-699 had the highest drug score with a value of 0.86. Furthermore, the toxicity properties evaluated were color-coded green, yellow or red. The properties displayed in red suggest a severe danger of unwanted consequences, yellow suggests mild toxicity, whereas the properties given in green imply drug conformity, compatibility, and safety in vivo. The toxicity results indicated that all the compounds showed no mutagenic or irritant risk. However, MolPort-047-733-181 was predicted to have high tumorigenic tendencies, while CSCS00058497692 with a medium irritant risk.
3.8. Molecular dynamics simulation
3.8.1. Root mean square deviation ligand and root mean square fluctuation
The stabilities of the c-alpha of the Pf 5-ALAS backbone without the ligand and in the protein-ligand complex as well as that of the CSMS00081585868 ligand in the complex were depicted in the RMSD plots obtained from the molecular dynamics (MD) simulation studies (Figures 6A–C). The RMSD of the two states of the c-alpha protein backbone fluctuated between 1.5 and 4.0 Å, with an average RMSD less that 3 Å. This confirmed that the conformation of the protein remained relatively stable throughout the MD simulation even with the binding of the ligand in its active site (Figure 6A). The RMSD of CSMS00081585868 fluctuated around 1.2 Å over the 10,020 frames after an initial increase from 0.0 Å, suggesting that the compound do not significantly change its orientation during the simulation (Figure 6B). The RMSF plot of the c-alpha protein backbone is shown in Figure 6C. A total of 630 residue positions were observed, while the highest RMSF values after simulation were observed at positions 118 and 119.
Figure 6. (A) RMSD plot of the c-alpha protein backbone (without ligand and in complex with ligand). (B) RMSD plot of CSMS00081585868 in the active site of Pf 5-ALAS. (C) RMSF plot of the c-alpha protein backbone.
3.8.2. Principal component analysis
The correlation of the motion of Pf 5-ALAS protein atoms are converted to a group of principal components which are uncorrelated, and this is captured in the form of PCA plots (Figure 7). The first PC (PC1) accounted for 53.5% of the cumulative variance, while PC2 and PC3 were responsible for 10.14 and 6.61%, respectively, as seen from the eigenvalue rank plot.
Figure 7. The principal component analysis (PCA) results for each data point (colored from blue to red according to time series) and the Eigenvalue rank plot.
3.8.3. Hydrogen bond analysis
H-bond formations are important in the stabilities of the complexes during MD simulation. The number of H-bonds observed between the atoms of CSMS00081585868 and the amino acid residues in the protein binding pockets during the simulation are shown in Figure 8. A total of 9 h-bonds was found during the MD simulation while the highest occupancy rate (19.93%) was observed for the interaction between CSMS00081585868 as the donor and ASP 430 as the acceptor. This corresponded with the observed h-hond interaction between the hydrogen of the NH on pyridine in CSMS00081585868 structural template and ASP 430 of Pf 5-ALAS at bond length of 2.20 Å from the docking studies (as seen in Figure 2B).
Figure 8. Hydrogen bonds formed between CSMS00081585868 in the active site of Pf 5-ALAS during the production run.
4. Discussion
Most of the available chemotherapeutic interventions for the treatment of malaria are designed against the erythrocytic stage, while only a few works have been done regarding the liver stage (41). This has been attributed to little knowledge about the biology of the malaria parasite, especially the deadly Plasmodium falciparum (42, 43). In the quest to get possible drug targets essential at the liver stage, Pf 5-ALAS model was developed using the ab initio modeling approach (44). The predicted model generated using the Robetta server (Figure 1) gave a better structure assessment result than the AlphaFold model (Table 2). Previous studies have also reported that in some cases, the Robetta model provided a better structural confidence score than the AlphaFold model (45). The docking scores from the virtual screening of downloaded compounds from nine different databases against the Pf 5-ALAS model showed that all the ten best hits possessed better binding affinity than the reported co-ligand, pyridoxal 5′-phosphate (6). The post-screening analysis carried out also helped to identify the binding interactions and bond lengths between the atoms of the best hits and the amino acid residues in the protein’s active site. CSMS00081585868 had the best ligand efficiency with a docking score of -9.9 kcal/mol and predicted Ki of 52.10 nM, engaging in seven hydrogen bonds with the target’s active site amino acid residues. The top ten hits had binding affinities in the order: CSMS00081585868 > MolPort-047-716-699 > CSCS00058497692 > MolPort-047-733-181 > ZINC5276997 > CSMS00083851644 > MolPort-039-017-789 > ZINC65533569 > MolPort-046-194-001 > MolPort-045-962-768 (Table 5). The docking scores can be attributed to the binding interactions observed between the atoms of the best hits and the binding pocket of the modeled target (24; Table 6).
The pharmacophore features of a compound go a long way to determining whether the compound will have good interactions with the target (46, 47). Heterocyclic compounds are important pharmacophores in drug design as over 90% of commercially-available drugs have them embedded in their structure, and they have a wide range of therapeutic applications (48). Hence, the four best hits from this study were examined to understand the structural components in their backbone responsible for the interactions observed (Figures 2–5). All the four best hits possessed heterocycles as major backbones in their structures, with the pyrrolidine template observed as a major linker in CSMS00081585868 (Figure 2) and MolPort-047-716-699 (Figure 3). The qualitative structural assessment of the best hit, CSMS00081585868, revealed its drug-like configuration, with two pyridine heterocyclic scaffolds bearing hydroxy and fluorine atoms linked by a pyrrolidine heterocyclic scaffold.
Drug attrition is one of the significant problems in medicinal research, as most proposed drug candidates at the discovery phase end up not being translated into marketable drugs (31). This has been attributed to poor pharmacokinetics and high toxic risk (49). Hence, in silico ADMET study was conducted on the best hits, and it was observed that all the compounds had good pharmacokinetic properties according to Lipinski’s rule of five (37). The drug scores of the ten best hits were in the order: MolPort-047-716-699 > CSMS00081585868 > CSMS00083851644 = ZINC65533 569 > MolPort-046-194-001 > MolPort-047-733-181 = ZINC 5276997 > CSMS00083851644 > MolPort-045-962-768 > CSCS00058497692. For the toxicity risk, MolPort-047-733-181 and CSCS00058497692 had high tumorigenic tendencies and medium irritant risk, respectively (Table 7). However, careful hit-to-lead optimization can be carried out on these compounds to make them non-toxic, as toxicity is usually a result of the presence of an unwanted functional group or pharmacophore present in the structure of the compound (50).
The stability of the best hit, CSMS00081585868, in the active site of Pf 5-ALAS and that of the c-alpha backbone of the protein were determined using MD simulation via the RMSD and RMSF analyses (Figure 6), PCA (Figure 7) and hydrogen bond analysis (Figure 8). Calculation of the RMSDs of the ligand and c-alpha protein backbone are important to identify the ligand’s possible modes of binding (51). The RMSD plot the CSMS00081585868 in the active site suggest that the ligand is stable with a steady fluctuation around 1.2 Å (52). The stabilities of both the ligand and protein during the MD simulation process showed that a stable complex was formed. A comparison of the c-alpha protein backbone without the ligand and in complex with the ligand showed that the RMSD observed from the complexed protein was due to the ligand binding. Inspection of the Pf 5-ALAS surface during simulation, using VMD, also showed that the major fluctuations especially those with RMSF above 1.0 Å were from the flexible loop regions (53). The variance proportion resulting from each principal component (PC) was also observed from the eigenvalue rank plot (Figure 7). The PCA study showed that the first three principal components accounted for 70.3% of the total variance while the hydrogen bond analyses (Figure 8) showed that the ligand maintained stable conformation in the active site of the protein during the simulation, suggesting the inhibitory potential of CSMS00081585868 against Pf 5-ALAS.
The limitation of this study is that little is still known about the biology of Plasmodium falciparum, especially the liver stage cycle of the parasite (42, 43). However, the advent of computer-aided drug design (CADD) has made it possible to predict the activities of several targets in the parasite and thereby design potential inhibitors against the targets. In conclusion, the essentiality of the heme biosynthetic pathway of Pf at the liver stage presents an opportunity to design prophylactic drugs to reduce the incidence and burden of malaria disease in endemic regions. In this study, following Pf 5-ALAS model design, target evaluation, macromolecular processing, library preparations of 2,621 compounds and processing, we employed structure-based virtual screening to identify potential lead compounds possessing higher inhibitory and drug-like properties to Pf 5-ALAS than the cofactor pyridoxal 5′-phosphate. We recommend that hit-to-lead (H2L) optimization should be employed, with the idea of the qualitative structural assessments, to improve the drug-like efficiency of the identified hits.
Data availability statement
The original contributions presented in this study are included in this article/supplementary material, further inquiries can be directed to the corresponding author.
Author contributions
EA, GO, RA, and IA spearheaded the project and collaborated with AV to design the work. GO, RA, and IA interpreted the results and wrote the first draft of the manuscript. All authors reviewed and approved the final draft.
Funding
This work was funded by the West Africa Sustainable Leadership and Innovation Training in Bioinformatics Research (WASLITBRe) through grant number 5U2RTW010679 and the Covenant Applied Informatics and Communication Africa Center of Excellence (2019-2025) World Bank.
Acknowledgments
The authors appreciate the Covenant University Bioinformatics Research (CUBRe) management for making available the necessary tools and infrastructure important for the successful execution of the project.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Footnotes
References
2. Menard D, Dondorp A. Antimalarial drug resistance: a threat to malaria elimination. Cold Spring Harb Perspect Med. (2017) 7:1–24. doi: 10.1101/cshperspect.a025619
3. Bazzani S, Hoppe A, Holzhütter H. Network-based assessment of the selectivity of metabolic drug targets in Plasmodium falciparum with respect to human liver metabolism. BMC Syst Biol. (2012) 6:118. doi: 10.1186/1752-0509-6-118
4. Fatumo S, Plaimas K, Adebiyi E, König R. Comparing metabolic network models based on genomic and automatically inferred enzyme information from Plasmodium and its human host to define drug targets in silico. Infect Genet Evol. (2011) 11:201–8. doi: 10.1016/j.meegid.2010.08.012
5. Nagaraj V, Sundaram B, Varadarajan N, Subramani P, Kalappa D, Ghosh S, et al. Malaria parasite-synthesized heme is essential in the mosquito and liver stages and complements host heme in the blood stages of infection. PLoS Pathog. (2013) 9:e1003522. doi: 10.1371/journal.ppat.1003522
6. Ikushiro H, Nagami A, Takai T, Sawai T, Shimeno Y, Hori H, et al. Heme-dependent Inactivation of 5-Aminolevulinate Synthase from Caulobacter crescentus. Sci Rep. (2018) 8:14228. doi: 10.1038/s41598-018-32591-z
7. Delves M, Plouffe D, Scheurer C, Meister S, Wittlin S, Winzeler E, et al. The activities of current antimalarial drugs on the life cycle stages of plasmodium: a comparative study with human and rodent parasites. PLoS Med. (2012) 9:e1001169. doi: 10.1371/journal.pmed.1001169
8. Goldberg D, Sigala P. Plasmodium heme biosynthesis: To be or not to be essential? PLoS Pathog. (2017) 13:e1006511. doi: 10.1371/journal.ppat.1006511
9. Rizopoulos Z, Matuschewski K, Haussig J. Distinct prominent roles for enzymes of Plasmodium berghei heme biosynthesis in sporozoite and liver stage maturation. Infect Immun. (2016) 84:3252–62. doi: 10.1128/IAI.00148-16
10. Carolino K, Winzeler E. The antimalarial resistome – finding new drug targets and their modes of action. Curr Opin Microbiol. (2020) 57:49–55. doi: 10.1016/j.mib.2020.06.004
11. de Sousa A, Combrinck J, Maepa K, Egan T. Virtual screening as a tool to discover new β-haematin inhibitors with activity against malaria parasites. Sci Rep. (2020) 10:3374. doi: 10.1038/s41598-020-60221-0
12. Anurak C, Kesara N. A systematic review: application of in silico models for antimalarial drug discovery. African J Pharm Pharmacol. (2018) 12:159–67. doi: 10.5897/AJPP2018.4904
13. Muratov E, Amaro R, Andrade C, Brown N, Ekins S, Fourches D, et al. A critical overview of computational approaches employed for COVID-19 drug discovery. Chem Soc Rev. (2021) 50:9121–51. doi: 10.1039/D0CS01065K
14. Ahmad S, Abbasi H, Shahid S, Gul S, Abbasi S. Molecular docking, simulation and MM-PBSA studies of nigella sativa compounds: a computational quest to identify potential natural antiviral for COVID-19 treatment. J Biomol Struct Dyn. (2021) 39:4225–33. doi: 10.1080/07391102.2020.1775129
15. Jumper J, Evans R, Pritzel A, Green T, Figurnov M, Ronneberger O, et al. Highly accurate protein structure prediction with AlphaFold. Nature. (2021) 596:583–9. doi: 10.1038/s41586-021-03819-2
16. Varadi M, Anyango S, Deshpande M, Nair S, Natassia C, Yordanova G, et al. AlphaFold protein structure database: massively expanding the structural coverage of protein-sequence space with high-accuracy models. Nucleic Acids Res. (2022) 50:D439–44. doi: 10.1093/nar/gkab1061
17. Chivian D, Kim D, Malmström L, Bradley P, Robertson T, Murphy P, et al. Automated prediction of CASP-5 structures using the Robetta server. Proteins Struct Funct Bioinforma. (2003) 53:524–33. doi: 10.1002/prot.10529
18. Santhoshkumar R, Yusuf A. In silico structural modeling and analysis of physicochemical properties of curcumin synthase (CURS1, CURS2, and CURS3) proteins of Curcuma longa. J Genet Eng Biotechnol. (2020) 18:24. doi: 10.1186/s43141-020-00041-x
19. Qi G, Tollefson M, Gogal R, Smith R, AlQuraishi M, Schnieders M. Protein structure prediction using a maximum likelihood formulation of a recurrent geometric network. bioRxiv [Preprint]. (2021) doi: 10.1101/2021.09.03.458873
20. Maiti S, Banerjee A. Epigallocatechin gallate and Theaflavin gallate interaction in SARS-CoV-2 spike-protein central channel with reference to the hydroxychloroquine interaction: Bioinformatics and molecular docking study. Drug Dev Res. (2021) 82:86–96. doi: 10.1002/ddr.21730
21. Tian W, Chen C, Lei X, Zhao J, Liang J. CASTp 3.0: computed atlas of surface topography of proteins. Nucleic Acids Res. (2018) 46:W363–7. doi: 10.1093/nar/gky473
22. Jendele L, Krivak R, Skoda P, Novotny M, Hoksza D. PrankWeb: a web server for ligand binding site prediction and visualization. Nucleic Acids Res. (2019) 47:W345–9. doi: 10.1093/nar/gkz424
23. Afolabi R, Chinedu S, Ajamma Y, Adam Y, Koenig R, Adebiyi E. Computational identification of Plasmodium falciparum RNA pseudouridylate synthase as a viable drug target, its physicochemical properties, 3D structure prediction and prediction of potential inhibitors. Infect Genet Evol. (2021) 97:105194. doi: 10.1016/j.meegid.2021.105194
24. Naqvi A, Mohammad T, Hasan G, Hassan M. Advancements in docking and molecular dynamics simulations towards ligand-receptor Interactions and Structure-function Relationships. Curr Top Med Chem. (2019) 18:1755–68. doi: 10.2174/1568026618666181025114157
25. Sunseri J, Koes D. Pharmit: interactive exploration of chemical space. Nucleic Acids Res. (2016) 44:442–8. doi: 10.1093/nar/gkw287
26. Prabitha P, Shanmugarajan D, Kumar T, Kumar B. Multi-conformational frame from molecular dynamics as a structure-based pharmacophore model for mapping, screening and identifying ligands against PPAR-γ: a new protocol to develop promising candidates. J Biomol Struct Dyn. (2022) 40:2663–73. doi: 10.1080/07391102.2020.1841677
27. Oduselu G, Ajani O, Ajamma Y, Adebiyi E. Structure-Based drug design in discovering target specific drugs against plasmodium falciparum adenylosuccinate lyase. Trop J Nat Prod Res. (2021) 5:739–43. doi: 10.26538/tjnpr/v5i4.23
28. O’Boyle N, Banck M, James C, Morley C, Vandermeersch T, Hutchison G. Open Babel: an open chemical toolbox. J Cheminform. (2011) 3:33. doi: 10.1186/1758-2946-3-33
29. Dallakyan S, Olson A. Small molecule library screening by docking with PyRx. Methods Mol Biol. (2015) 1263:1–11.
30. Sravani M, Kumaran A, Dhamdhere AT, Senthil Kumar N. Computational molecular docking analysis and visualisation of anthocyanins for anticancer activity. Int J Res Appl Sci Biotechnol. (2021) 8:154–61. doi: 10.31033/ijrasb.8.1.18
31. Agamah F, Mazandu G, Hassan R, Bope C, Thomford N, Ghansah A, et al. Computational/in silico methods in drug target and lead prediction. Brief Bioinform. (2020) 21:1663–75. doi: 10.1093/bib/bbz103
32. Oduselu G, Ajani O, Ajamma Y, Brors B, Adebiyi E. Homology Modelling and Molecular Docking Studies of Selected Substituted Benzo[d]imidazol-1-yl)methyl)benzimidamide Scaffolds on Plasmodium falciparum Adenylosuccinate Lyase Receptor. Bioinform Biol Insights. (2019) 13:1177932219865533. doi: 10.1177/1177932219865533
33. Phillips J, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, et al. Scalable molecular dynamics with NAMD. J Comput Chem. (2005) 26:1781–802. doi: 10.1002/jcc.20289
34. Humphrey W, Dalke A, Schulten K. VMD: Visual Molecular Dynamics. J Mol Graph. (1996) 7855:33–8. doi: 10.1016/0263-7855(96)00018-5
35. Grant B, Rodrigues A, Elsawy K, Mccammon J, Caves L. Bio3d: an R package for the comparative analysis of protein structures. Bioinformatics. (2006) 22:2695–6. doi: 10.1093/bioinformatics/btl461
36. Thillainayagam M, Malathi K, Ramaiah S. In-Silico molecular docking and simulation studies on novel chalcone and flavone hybrid derivatives with 1, 2, 3-triazole linkage as vital inhibitors of Plasmodium falciparum dihydroorotate dehydrogenase. J Biomol Struct Dyn. (2017) 36:3993–4009. doi: 10.1080/07391102.2017.1404935
37. Lipinski C, Lombardo F, Dominy B, Feeney P. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Deliv Rev. (1997) 23:3–25. doi: 10.1016/S0169-409X(96)00423-1
38. Husain A, Ahmad A, Khan S, Asif M, Bhutani R, Al-Abbasi F. Synthesis, molecular properties, toxicity and biological evaluation of some new substituted imidazolidine derivatives in search of potent anti-inflammatory agents. Saudi Pharm J. (2016) 24:104–14. doi: 10.1016/j.jsps.2015.02.008
39. Gogoi P, Shakya A, Ghosh S, Gogoi N, Gahtori P, Singh N, et al. In silico study, synthesis, and evaluation of the antimalarial activity of hybrid dimethoxy pyrazole 1,3,5-triazine derivatives. J Biochem Mol Toxicol. (2021) 35:1–13. doi: 10.1002/jbt.22682
40. Behrouz S, Navid M, Rad S, Taghavi B, Mohammad S, Behrouz M, et al. Design, synthesis, and in silico studies of novel eugenyloxy propanol azole derivatives having potent antinociceptive activity and evaluation of their β -adrenoceptor blocking property. Mol Divers. (2018) 23:147–64. doi: 10.1007/s11030-018-9867-7
41. Goswami D, Minkah N, Kappe S. Malaria parasite liver stages. J Hepatol. (2022) 76:735–7. doi: 10.1016/j.jhep.2021.05.034
42. Chugh A, Kumar A, Verma A, Kumar S, Kumar P. A review of antimalarial activity of two or three nitrogen atoms containing heterocyclic compounds. Med Chem Res. (2020) 29:1723–50. doi: 10.1007/s00044-020-02604-6
43. Birnbaum J, Flemming S, Reichard N, Soares A, Mesén-Ramírez P, Jonscher E, et al. A genetic system to study Plasmodium falciparum protein function. Nat Methods. (2017) 14:450–6. doi: 10.1038/nmeth.4223
44. Kaushik A, Kumar S, Wei D, Sahi S. Structure based virtual screening studies to identify novel potential compounds for GPR142 and their relative dynamic analysis for study of type 2 Diabetes. Front Chem. (2018) 6:23. doi: 10.3389/fchem.2018.00023
45. Goh N, Mohamad Razif M, Yap Y, Ng C, Fung S. In silico analysis and characterization of medicinal mushroom cystathionine beta-synthase as an angiotensin converting enzyme (ACE) inhibitory protein. Comput Biol Chem. (2022) 96:107620. doi: 10.1016/j.compbiolchem.2021.107620
46. Marinescu M. Synthesis of Antimicrobial Benzimidazole – Pyrazole Compounds and Their Biological Activities. Antibiotics. (2021) 10:1–29. doi: 10.3390/antibiotics10081002
47. Sharma V, Gupta M, Kumar P, Sharma AA. Comprehensive Review on Fused Heterocyclic as DNA Intercalators: Promising Anticancer Agents. Curr Pharm Des. (2021) 27:15–42. doi: 10.2174/1381612826666201118113311
48. Ajani O, Aderohunmu DV, Olorunshola S, Ikpo C, Olanrewaju I. Facile Synthesis, Characterization and Antimicrobial Activity of 2-Alkanamino Benzimidazole Derivatives. Orient J Chem. (2016) 32:109–20. doi: 10.13005/ojc/320111
49. Ayati A, Falahati M, Irannejad H, Emami S. Synthesis, in vitro antifungal evaluation and in silico study of 3-azolyl-4-chromanone phenylhydrazones. Daru. (2012) 20:1–7. doi: 10.1186/2008-2231-20-46
50. Wu J, Jiang R, Lin W, Ouyang G. Effect of salinity and humic acid on the aggregation and toxicity of polystyrene nanoplastics with different functional groups and charges. Environ Pollut. (2019) 245:836–43. doi: 10.1016/j.envpol.2018.11.055
51. Bray S, Senapathi T, Barnett C, Grüning B. Intuitive, reproducible high-throughput hmolecular dynamics in Galaxy: a tutorial. J Cheminform. (2020) 12:1–13. doi: 10.1186/s13321-020-00451-6
52. Aksimentiev A, Arkhipov A, Birnbaum R, Brunner R, Cohen J, Dhaliwal B, et al. Using VMD. Chicago, IL: University of Illinois (2019).
53. Johnson T, Adegboyega A, Ojo O, Yusuf A, Iwaloye O, Ugwah-Oguejiofor C, et al. A computational approach to elucidate the interactions of chemicals from Artemisia annua targeted toward SARS-CoV-2 main protease inhibition for COVID-19 treatment. Front Med. (2022) 9:907583. doi: 10.3389/fmed.2022.907583
Keywords: 5-ALAS, drug design, heterocycles, malaria, pharmacophore modeling, structure assessment
Citation: Oduselu GO, Afolabi R, Ademuwagun I, Vaughan A and Adebiyi E (2023) Structure-based pharmacophore modeling, virtual screening, and molecular dynamics simulation studies for identification of Plasmodium falciparum 5-aminolevulinate synthase inhibitors. Front. Med. 9:1022429. doi: 10.3389/fmed.2022.1022429
Received: 18 August 2022; Accepted: 23 December 2022;
Published: 12 January 2023.
Edited by:
Duoquan Wang, Chinese Center For Disease Control and Prevention, ChinaReviewed by:
Ana Carolina Rennó Sodero, Federal University of Rio de Janeiro, BrazilMithun Rudrapal, Rasiklal M. Dhariwal Institute of Pharmaceutical Education and Research, India
Copyright © 2023 Oduselu, Afolabi, Ademuwagun, Vaughan and Adebiyi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ezekiel Adebiyi, ZXpla2llbC5hZGViaXlpQGNvdmVuYW50dW5pdmVyc2l0eS5lZHUubmc=
†These authors have contributed equally to this work and share first authorship