Molecular Docking, Pharmacokinetics and Molecular Dynamics Simulation Studies of Some Bioactive Compounds Isolated from Entandrophragma congoënse for Antiplasmodial Activity

As a follow-up to earlier reported works on the phytochemical study of some isolated bioactive compounds from the root and bark of Entandrophragma congoënse as potent anti-plasmodium drugs (Happi et. al. 2005), some of the isolated compounds were tested in vitro for antiplasmodial and cytotoxicity but no insight was given into the binding


INTRODUCTION
Malaria is an infectious disease caused by the transmission of one of the five Plasmodium species including P. falciparum, P. ovale, P. malariae, P. knowlesi and P. vivax, through the bite of mosquitoes female Anopheles from an infected person to another one [1 ]. The disease is still classified among the most harmful that affects the population in the tropical and subtropical regions of the world [2]. A recent report by World Health Organization (WHO) on malaria indicated that in 2020, there were an estimated 241 million malaria cases and 627,000 deaths in 85 malaria-endemic countries with an increased case incidence due to the disruption to services during the COVID-19 pandemic. More particularly, six African countries including Nigeria (27%), the Democratic Republic of the Congo (12%), Uganda (5%), Mozambique (4%), Angola (3.4%) and Burkina Faso (3.4%) accounted for about 55% of all malaria cases and deaths globally in 2020 [3 ]. Several strategies developed in the last decades to tackle the progress of the disease include the large distribution of long-lasting insecticide-treated bed nets and the development of new insecticides as well as new potent drugs [2]. Numerous reported cases of the resistance of P. falciparum to many antimalarial drugs like quinine and artemisinin [4,5] necessitate a continuous search and development of new chemotherapeutic compounds to address the current situation of drug resistance [6]. Traditional medicine is an important part of health care in several countries of the world including the use of medicinal plants for the treatment of several diseases like malaria or more recently COVID-19 [7,8 ]. Relentless efforts in the search for new potent antiplasmodial metabolites from medicinal plants led to the isolation of numerous specialized metabolites from Cameroonian medicinal plants including a series of protolimonoids displaying a strong antiplasmodial activity against P. falciparum NF54. [1,9,10]. In order to gain further insights into their inhibitory mechanism of action against malaria parasites, their drug druglikeliness, toxicity, and other pharmacokinetic properties, the most active compounds isolated from Entandrophragma congoënse are computationally studied using molecular docking, ADMET studies and Molecular dynamics simulation techniques.

Ligand Preparation
In a quest to study the inhibitory effects of the selected ligands as anti-plasmodium potent drugs, all the twenty-one compounds (including two standards) were sketched and cleaned using the ACD ChemSketch ® 2.0 software [11]. Merck Molecular Force Field (MMFF) was adopted for conformation search using the Spartan 14 ® software 1.1.4 version [12] in order to obtain the most stable conformers. The Density Functional Theory method (DFT) at B3LYP and 6−31+G (d) basis set was used for geometry optimization of the most stable conformers.

Preparation of Target Receptor
The appropriate antiplasmodial receptor with the PDB code 5TBO was obtained from the Protein Data Bank (www.rcsb.org) web server [13] for this study. Water molecules around the crystal structure and other associated inhibitors were removed from the downloaded crystal complex to avoid undesired interferences and unwanted molecular interactions. Binding pockets of the complex were obtained from the associated literature in the RCSB database [14] and were confirmed computationally using the Computed Atlas for Surface Topography of proteins (CASTp) web server [15]. The molecular interactions were obtained using the BIOVIA Discovery Studio (2021) software [16 ].

Molecular Docking and Molecular Dynamics Simulation Studies
The molecular docking was done using the Pyrx/Vina tool [17]. and AutoDock Tools ® 1.5.6 version [18] and the ligands and macromolecules (in PDB format) were saved as (.pdbqt) file in preparation for docking. The grid box (X = 40.552, Y = 18.551 and Z = 84.308) was set in response to the determined active sites. Adjustments done to improve the affinity of the binding sites of the crystal structures includes removing any extra water molecules, adding polar hydrogens, merging no-polar hydrogens, and adding Kollman charges. The active sites of the studied antiplasmodial receptor (PDB code 5TBO) are VAL532, ILE237, LEU531, HIE185, TYR528, ASN274, ARG265. Molecular Dynamics Simulation studies were carried out by employing NAMD 2.0 version [19] with the minimum energy configurations of the compounds using CharmmFF forcefield. For ligands parameterization and topology, CHARMM GUI web server (www.charmmgui.org/) [20 ] was used. Solvation of the system was done via Simple Point Charge model (SPC/E) water, in a cubic box leaving 5.0 nm space around the solute. Counter ions in the form of Na + at a concentration, 0.15 M ions were added to neutralize all the systems.
Energy minimization of the system was carried out to reach a maximum force. The system was then equilibrated for 1000 ps at 300 K using an Isochoric-Isothermal or (NVT) ensemble followed by Isothermal-Isobaric or (NPT) ensemble for another 1000 ps of equilibration at 300 K. For both NVT and NPT equilibrations, the electrostatic and van der Waals interaction cutoffs were fixed at 1.0 nm. These equilibrated ensembles were then subjected to Molecular Dynamics simulation for 70 ns with the same electrostatic and van der Waals cutoffs.
The trajectories and structures of the complexes were observed and visualized by VMD (Visual Molecular Dynamics) [21]. For graphical analysis, the Xmgrace ® tool [22 ] was utilized to plot the graphs of Root Mean Square Deviation and Fluctuation (i.e. RMSD and RMSF), hydrogen bonding and Radius of Gyration (RoG).

ADMET and Drug-likeness Studies
To predict the toxicity level of all the twenty-one compounds, in-silico toxicity studies were performed using the admetSAR 2.0 web server [23 ] and the SwissADME ® tool [24 ], where the bio-availability radar was obtained. From the server, the drug-induced hERG toxicity, carcinogenicity, human oral bio-availability, AMES mutagenesis, acute oral toxicity, watersolubility, bio-degradation and other parameters were computed. Important conditions for a drug ADMET profile translates that a drug molecule should result to a good Human Intestinal Absorption (HIA), Solubility (Log S) range between -1 and -7, should not be an inhibitor of cytochrome P450 enzymes, and should be non-Ames toxic. Others include non-carcinogenicity, non-inhibition of hERG, and no or low level of toxicity [23,25 ].      Table 1 shows the binding energies and inhibition constants of the studied 21 compounds (2 standards inclusive) with 5TBO receptor. From the results, only compounds L1, L2, L4 and L15 fall within the recommended range against 5TBO receptor, and their molecular interactions with the binding site of the receptor are shown in Table 2. Hence, they can be considered as potent hits against Plasmodium. Furthermore, none of the standards (i.e. Azadirachtin ® and Chloroquine ® ) is qualified, as both fail to fall between recommended range values with 5TBO receptor. However, the study proceeded with the inclusion of compound L21 (Chloroquine ® ) as a standard together with the 4 qualified compounds in further studies, as potential anti-plasmodium targets.

Molecular Docking and Molecular Dynamics simulation studies
To further confirm the docking results, Table 3 shows the interactions formed by each qualified compound at the binding site of the 5TBO receptor. The high binding affinity score of compound L15 is justified by the short length of the conventional hydrogen bonding and other intermolecular interactions. The low binding affinity or rather, the high energy score of compound L21 (standard) is justifiable considering the high number of other weak interactions such as van der Waals, ℼ→ℼ stacked, ℼ→alkyl, ℼ→ and ntera t on instead of the stronger conventional hydrogen bond.
To gain more insight into structural mobility and atom flu tuat on from the y tem ' n t al coordinates over simulation trajectories, molecular dynamics studies were conducted. Fig. 1 shows the Root Mean Square Deviation (RMSD) of the 4 compounds with the standard. A total of 70 ns MD simulations were performed in order to evaluate the mulat on' on ten y, root mean quare deviation of backbone L1, L2, L4, L15, L21 and our target protein were analyzed and plotted. In Fig. 1, the RMSD curves of the protein backbone atoms show that the complexes had equilibrated and were stable after 70 ns. The RMSD average values from the plot for L1, L2, L4, L15, L21 and our target protein were 0.38, 0.40, 0.28, 0.22, 0.42 and 0.50 Å, respectively. The sequence of stability strength is decreasing in L15, L4, L1, L2, L21 and protein target which means that the compound L15 is the most stable in the study. So the RMSD plot shows the reliability of the MD simulation equilibrium and the stability of complexes after 70 ns.
The Root Mean Square Fluctuations (RMSF) of the alpha carbon atoms were computed to provide explanations and give a better understanding of the flexibility and structural fluctuations of the different regions of amino acids in L1, L2, L4, L15, L21 and our target protein complexes as shown in Fig. 2. From the results, it is clear that there are fluctuations in the residues of all systems within the range of 40 and 50; and from 105 to 215 amino acid residues. In these two regions, the fluctuations of the complexes are present. In the two fluctuated regions, we observed that compound L21 (red) has high fluctuations in the residues among all compound complexes while other complexes are less fluctuated.
The radius of gyration reveals the proteincomplex structural compactness. The higher radius of gyration value corresponds to a less tightly packed system and less stable conformation [26 ]. The average radius of gyration values of all systems such as L1, L2, L4, L15, and L21 are between 19.8 to 20Å as shown in Fig. 3. This implies that all the systems have approximately the same radius of gyration value which means that all the systems are tightly packed and stable. However, since, L1 has the lowest value of 19.88 among all, it is convenient to say that L1 is tightly packed and stable in terms of the radius of gyration analysis.

Hydrogen Bonds Analysis
Hydrogen bonds are considered as strong dipole-dipole interactions. Hydrogen bonds and their number are more important for the proteinligand complex. The stability of the ligand in the active binding cavity of protein is measured in terms of the average number of hydrogen bonds and the length of the interactions. The shorter the length, the shorter will be the interaction and the more stable the hydrogen bond will be [27]. Fig.  4 shows that all the qualified compounds retained hydrogen bonds over 200 after the minimization run of 50 ns and simulation run of 70 ns whereas, the standard, L21, fell short below 200 through the simulation period. It was revealed in Fig. 4 that the complex L15 (shown in black) has a high number of hydrogen bonds, at 2.5 ns simulation time the number of hydrogen bonds increases to 148 while in the rest of the simulation time, it is around 135. All the system except L15 has less number of hydrogen bonds throughout all simulation time. This result also confirms that L15 is the ligand that is holding the strongest to the binding site of the protein, hence revalidating the docking results.

ADMET and Drug-likeness Studies
Toxicity studies are needed to determine the possible adverse effects of a drug in humans, animals, plants and the environment. Table 4 shows the ADMET properties of the 4 qualified compounds and the standard (Chloroquine ® , L21). The results show that L21 is an inhibitor of (3) CYP450 which is significantly against the condition for a drug as reported by (Tsaioun and Kates 2011) . All selected compounds also possess a good acute oral toxicity rating of III. This can however be further optimized to IV during experimental considerations. All (4) selected compounds are non-carcinogenic, non-hERG toxic, have high gastrointestinal absorption and are likely to be used as a drug with L15 appearing with the highest probability (oral bioavailability of 0.85). The bioavailability radar further helps to show the physicochemical rating for oral bioavailability. The charts in Fig. 5 show that all the selected compounds, likewise, the standard, L21, are slightly within the coloured acceptable range with L15 diverting insignificantly due to its insolubility and lipophilicity which can be optimized to an acceptable range in further experimental analysis.

CONCLUSION
In this study, attempts are made to computationally gain insight into the binding affinities of the 19 isolated bioactive compounds from the root and bark of Entandrophragma congoënse as potent anti-plasmodium drugs using molecular docking and molecular dynamics simulation studies, studies, as well as ADMET and drug-likeness. Results of molecular docking studies revealed that only 4 (coded L1, L2, L4 and L15) of the 19 compounds are qualified based on the values of their binding energies and molecular interactions. ADMET and druglikeness studies show that all the 4 qualified compounds are non-tox and obey L p n k ' rules of drug-like compounds. Molecular dynamics studies confirmed the stabilities of the complexes of L1, L2, L4 and L15 (in comparison with the Chloroquine ® ). The overall results reveal that compound L15, when isolated, can alone, or together with other qualified compounds is/are better drug candidate(s) for the treatment of malaria compared to the FDA-approved drug (Chloroquine ® , L21).