Molecular Docking Investigation, Pharmacokinetic Analysis, and Molecular Dynamic Simulation of Some Benzoxaborole-Benzimidazole Hybrids: An Approach to Identifying Superior Onchocerca Inhibitors

Onchocerciasis is one of the major neglected tropical diseases caused by the filarial worm ( Onchocerca volvulus ), affecting an estimated population of about 37 million people living predominantly in tropical Africa. The major treatment approach has been based on the use of Ivermectin, which kills the microfilariae or the less effective Doxycycline targeting Wolbachia , endosymbiont of filarial nematodes. Flubendazole (FBZ) has proved effective in treating adult worms but with threatening adverse effects. Against this backdrop, therefore, a combined molecular docking study and pharmacokinetic screening were conducted on a series of benzimidazole-benzoxaborole hybrids to find more potent analogs with attributes that address the limitations of existing therapies. All the nineteen analogs were found to possess better docking scores than the reference drug (FBZ, Moldock scores = -120.466 and - 125.359). The results of pharmacokinetic testing suggest that four molecules ( 14 , 16 , 19 , and 20 ) are orally bioavailable and showed better ADMET properties than FBZ. These molecules and FBZ showed good binding interactions with the receptors’ active sites. Also, the molecular dynamic simulation performed on the docked complexes of 20 and FBZ confirmed the rigidity and stability of their interactions. Based on the results of this study, the selected molecules (especially 20 ) could be considered superior drug candidates for the treatment of Onchocerciasis.


INTRODUCTION Abstract
Onchocerciasis is one of the major neglected tropical diseases caused by the filarial worm (Onchocerca volvulus), affecting an estimated population of about 37 million people living predominantly in tropical Africa. The major treatment approach has been based on the use of Ivermectin, which kills the microfilariae or the less effective Doxycycline targeting Wolbachia, endosymbiont of filarial nematodes. Flubendazole (FBZ) has proved effective in treating adult worms but with threatening adverse effects. Against this backdrop, therefore, a combined molecular docking study and pharmacokinetic screening were conducted on a series of benzimidazole-benzoxaborole hybrids to find more potent analogs with attributes that address the limitations of existing therapies. All the nineteen analogs were found to possess better docking scores than the reference drug (FBZ, Moldock scores = -120.466 and -125.359). The results of pharmacokinetic testing suggest that four molecules (14,16,19, and 20) are orally bioavailable and showed better ADMET properties than FBZ. These molecules and FBZ showed good binding interactions with the receptors' active sites. Also, the molecular dynamic simulation performed on the docked complexes of 20 and FBZ confirmed the rigidity and stability of their interactions. Based on the results of this study, the selected molecules (especially 20) could be considered superior drug candidates for the treatment of Onchocerciasis.

Data sourcing
Akama et al. 9 and Carter et al. 10 reported the synthesis of a series of amide-linked and ketone-linked analogs of benzoxaborole-benzimidazole hybrid, respectively, as well as FBZ, as part of the anti-Onchocerca drug discovery effort. Their biological activities were tested against O. volvulus L3 larval molting assay and reported in micromolar (µM). Consequently, a dataset of 20 compounds comprising FBZ and benzoxaborole-benzimidazole analogs with relatively better half-maximal inhibitory concentration (IC50) values was obtained from their report and used for this theoretical study. The molecular structures, bioactivities (IC50), and pIC50 obtained as a logarithmic function of IC50 for these compounds were reported in Table I.

Structural optimization
The molecular structures of all the compounds in Table I were drawn using the ChemDraw Ultra, saved as MDL molfile format and fed separately onto the Spartan '14 Graphical User Interface. Energy minimization was performed on the imported molecules and then saved in Spartan file format. The resulting structures were then subjected to full-scale optimization first by using Molecular Mechanics Force Field (MMFF) and after that, Density Functional Theory (DFT) with Becke's three-parameter read-Yang-Parr hybrid (B3LYP) option and utilizing the 6-31G** basis set. The optimized structures were saved as PDB formats for subsequent use in molecular docking studies 13,24 .

Docking protocol
The crystal structures of two O. volvulus target proteins, prostaglandin D synthase (PDB: 2HNL) and pi-class glutathione Stransferase (PDB: 1TU7), were retrieved from the RCSB Protein Data Bank in the PDB file format and then modified separately using the Molegro Virtual Docker (MVD) by eliminating water molecules, cofactors, and co-crystallized ligands contained within the protein structures 25 . Both proteins comprise similar chains (A and B), while A was utilized. The software allows for the repair (rebuild) of all affected residues. The receptor's binding cavities were defined, and those which have the largest volume and surface areas (volume: 38.912 and 33.28; surface: 139.52 and 144.64) for 2HNL and 1TU7 respectively, were adopted for the docking. All ligands were imported in PDB file format and prepared appropriately. The simulation was performed using the parameter settings in Table II. The binding scores were then recorded, while the predicted ligand-protein interaction profiles were visualized using the Biovia Discovery Studio Visualizer. A similar method was earlier reported elsewhere 14, 26 .

Prediction of pharmacokinetic properties
Drug-likeness and ADMET properties prediction constitute a vital stage in drug discovery's early phase because only molecules with good drug-likeness properties and excellent ADMET profiles advance into the pre-clinical research phase 13 .
Here, FBZ and the selected compounds (6, 14, 15, 16, 18, 19, and 20) with the best binding scores were investigated for their pharmacokinetic properties using three online web servers: SwissADME, pkCSM, and admetSAR for drug-likeness, ADME and toxicity profiling respectively. The choices of molecules for oral bioavailability have been guided by several rules such as Lipinski's 'rule of 5' (RO5), Veber rule, Ghose rule, Egan, Muegge, and others 27 . Lipinski's RO5 is a widely used criterion for oral bioavailability. As such, these compounds would be assessed for oral bioavailability using the RO5 criteria 15 .

Molecular dynamics simulation
Molecular dynamics (MD) simulations of the complexes of compound 20 with both proteins (2HNL and 1TU7) were performed separately using the combined approach of Chemistry at Harvard Macromolecular Mechanics (CHARMM) force field, Nano-scale Molecular Dynamics (NAMD), and Visual Molecular Dynamics (VMD). The CHARMM-GUI, an established web-based platform that utilizes the CHARMM force field, was used to generate the input files for the simulation by NAMD 23 . The periodic boundary condition was utilized while fitting the system into a cubic water box for solvation. The protein was solvated and neutralized explicitly in an aqueous solution of 0.1 M KCl salt 16 . Energy minimization was performed to stabilize the complex structure and ensure steric clashes will not result. The resulting system of ions and solvent was then equilibrated to stabilize the system at a temperature chosen for the simulation (310 K) at a constant number of particles, volume, and temperature (NVT ensemble) and to stabilize the pressure by keeping the number of particles, pressure, and temperature (NPT ensemble) constant using 100 ps time frame 28 . MD was then performed on the resulting system for 10 ns (5,000,000 steps), while the results were visualized using VMD and the Biovia discovery studio. A similar procedure was described elsewhere 16 . Additionally, MolAICal software was used to compute the ligand-binding affinity by Molecular Mechanics Generalized Born Surface Area (MM/GBSA) method based on the resulting MD log files obtained with NAMD 29 . MM/GBSA is estimated using Equations 1 to 3, where ∆EMM is the gas phase MM energy; -T∆S represents the conformational entropy; ∆EMM is the sum of ∆Eele, van der Waals energy ∆Evdw and ∆Einternal of bond, angle, and dihedral energies; ∆Gsol is the solvation free energy equal to the sum of the non-electrostatic solvation component ∆GSA and electrostatic solvation energy ∆GGB.

Virtual docking screening
The results (Moldock scores) of the docking simulation conducted between the two receptors of O. volvulus and the various benzimidazole-benzoxaborole hybrids, as well as the reference drug (flubendazole), were reported in Table III, while Figure 1 shows the 3D representation of compound 20 in the active sites of both target proteins. The binding score, which indicates the affinity between a ligand and a receptor, is often used to screen an extensive library of compounds to find more active molecules interacting strongly with the receptor of interest. The binding affinities (Moldock scores) available in Table  III  More so, compound 20 is promising as it binds relatively very strongly with both receptors (characteristic of multi-target drug molecules) while having the highest reported biological activity (pIC50) of 8.3979, the same as that of FBZ. Therefore, the virtual docking screening has proved effective as the most active analogs were identified and extracted for further evaluation.

Evaluation of pharmacokinetic properties
Drug-likeness analysis and ADMET studies were conducted on FBZ and the seven selected analogs (6, 14, 15, 16, 18, 19, and 20) to ascertain their oral bioavailability, toxicity, and safety profiles. The results of both investigations were presented in Tables IV and V, respectively, with Figures 2 and 3 showing the bioavailability radar, while Figure 4 showed their Boiled Egg's representation. Lipinski's approach to ascertaining the oral bioavailability of compounds has been widely applied in the discovery of new drug molecules 13 . It asserts that a drug molecule may likely not be bioavailable orally when it has an HBD greater than 5, HBA > 10, MW > 500, and MLOGP > 4.15 or WLOGP > 5 15 . Whenever a molecule passes at least three of the four provisions of the RO5, it is said to comply with Lipinski's rule for oral bioavailability 14 . Table IV shows that FBZ and all the tested benzimidazole-benzoxaborole analogs passed the drug-likeness test (Lipinski RO5) by showing no violation. The reported Topological Polar Surface Area (TPSA) values for the molecules were below the threshold value of 140 Å 2 , beyond which a molecule may exhibit poor Human Intestinal Absorption (HIA). Also, the synthetic accessibility (SA) scores of these compounds were less than 5.00 (accessible portion on a scale of 1 to 10), suggesting easy laboratory synthesis of these molecules. The predicted values of the estimated water solubility (Log S) range of -2 > Log S >-4 indicate that these molecules are aqueous soluble. None of the compounds were also predicted as Pan Assay Interference compounds (PAINS). The estimated ADMET properties reported in Table V showed good Human Intestinal Absorption (HIA) (greater than 65%) for all tested compounds. Skin permeability is a key factor in transdermal drug delivery development, with skin permeation constant LogKp >-2.50 indicating poor skin permeability. As a result, the various compounds showed LogKp values <-2.50, connoting good skin permeability. Drug molecule penetration through the Blood-Brain Barrier (BBB) and Central Nervous System (CNS) comes with specific criteria, which specify that for a drug molecule to penetrate the BBB and CNS readily, the logarithmic ratio of brain-to-plasma drug concentration (logBB) must be >0.3 and the blood-brain permeability-surface area product (logPS) be >-2 respectively. Consequently, none of these molecules were predicted to penetrate the BBB and CNS readily. Furthermore, some groups of cytochromes P450 enzymes are important in the body to facilitate drug metabolism and help in their excretion. The two major isoforms enhancing drug metabolism, CYP34A and CYP2D6, were tested. However, the tested molecules were not substrates and inhibitors of both enzymes. Also, FBZ only is a non-substrate of P-glycoprotein, an enzyme that acts as a biological barrier by extruding toxins and xenobiotics, including drugs, out of cells. This means that FBZ may not be effluated out of the target cells by this enzyme when taken into the human system. However, like FBZ, only 6, 14, and 20 showed inhibition of this enzyme. The drug's total clearance determines the extent of drug removal from the body. The range of total clearance values for all the tested molecules is good. The renal Organic Cation Transporter 2 (OCT2) is critical in drug and endogenous compounds' disposition and renal clearance. However, the activities of OCT2 may cause unwanted drug-drug interactions, which makes it necessary to ascertain whether a drug molecule is an OCT2 substrate or inhibitor. As seen in Table V, all the tested compounds are non-OCT2 substrates. Additionally, some toxicity indices were predicted to ascertain the safety profiles of these molecules. The Ames test is widely applied to ascertain a compound's mutagenic potential. A positive test indicates that the compound is mutagenic. Only compounds 6, 15, and 18 showed positive Ames toxicity, showing a mutagenicity risk. Also available in Table V is the Maximum Recommended Tolerated Dose (MRTD) predicted for the various molecules. MRTD value of ≤ 0.477 log (mg/kg/day) is considered low, while a value > 0.477 log (mg/kg/day) is considered high. The inhibition of the human ether-a-go-go gene (hERG) is responsible for the acquired long QT syndrome, resulting in heartbeat irregularity issues. Only FBZ was implicated in this regard. Furthermore, all the tested compounds showed positive hepatotoxicity but negative skin sensitization, eye irritation, and carcinogenicity. Additionally, only FBZ and compound 6 are nephrotoxic, indicating the risk of rapid deterioration in kidney function. Based on the toxicity analysis, compounds 14, 16, 19, and 20 exhibited relatively safer toxicity profiles. Therefore, further evaluations were limited to these molecules and FBZ. The bioavailability radar displayed in Figures 2 and 3 was to provide a rapid appraisal of the molecules' drug-likeness. The radar considers six physicochemical properties: size, polarity, lipophilicity, solubility, flexibility, and saturation. The radar's pink area represents the ideal or suitable physicochemical space in which the molecule falls completely to be classified as drug-like. This means lipophilicity (-0.7 ≤ XLOGP3 ≤ +5.0), size (150 ≤ MW ≤ 500), polarity (20 ≤ TPSA ≤ 130), solubility (Log S ≤6), saturation (fraction of carbons in the sp3 hybridization ≥0.25), and flexibility (No. of rotatable bonds ≤9) 22 . The selected molecules (14, 16, 19, and 20) were said to be orally bioavailable, as seen from their bioavailability radar in Figure 2. FBZ, on the other hand, showed high unsaturation due to its low fraction of carbon in the sp3 hybridization and, therefore, is not orally bioavailable. This conforms with reports available in the literature regarding the need to improve the oral bioavailability of FBZ 8 . Figure 4 shows the boiled egg representation of FBZ and the four molecules of benzimidazole-benzoxaborole analogs. All the molecules were located in the boiled egg's white, indicating that they were predicted to be passively absorbed by the gastrointestinal tract. Also, only the reference drug appeared in the red dot, showing that it was predicted not to be effluated from the CNS by the pglycoprotein. The four compounds' overall drug-likeness and ADMET properties (14, 16, 19, and 20) showed better pharmacokinetic profiles than the reference drug, especially in bioavailability and toxicity.

Pharmacological interactions
The pharmacological interactions between the receptors' amino acid residues and the selected compounds (14, 16, 19, and 20), as well as the reference drug (FBZ), were summarized in Table VI, while the 2D and 3D views of the binding interactions as adapted from the Discovery Studio Visualizer were shown in Figures 5 to 14. This was to provide insight into the mode of binding of these ligands with the active sites of the various target proteins. The various compounds were said to interact adequately with the respective target receptors, as shown by the presence of hydrogen bonding (H-bond), electrostatic interactions, and hydrophobic interactions (Table VI). H-bonding and hydrophobic interactions were present in all the interactions involving the two target proteins, while electrostatic interactions were more visible in the interactions involving 1TU7 than in those of 2HNL. Incorporating the benzoxaborole group into the benzimidazole core significantly improved the binding interactions of the selected analogs with the receptors, as seen in the higher binding affinities associated with these molecules. Because compound 20 interacted very strongly with both receptors (Moldock score= -154.024 for 2HNL and -157.696 for 1TU7) and also had the highest pIC50 value of 8.3979, its binding interactions were discussed and compared with those of the reference drug (FBZ). A total of two conventional H-bonds, two carbon H-bonds, two π-anion electrostatic interactions, and up to eight hydrophobic interactions were formed between compound 20 and 1TU7 (Figure 12). On the other hand, the binding interaction profile of FBZ with 1TU7 involved four conventional H-bond interactions, two π-cation and two π-anion electrostatic interactions, and up to five hydrophobic interactions (Figure 14). Unlike in 20_1TU7, no unfavorable steric clash was visible. The conventional H-bonds were formed with ARG-195 at a distance of 3.10 Å, LYS-183 (2.33Å), and VAL-200 at distances of 2.20 and 2.89 Å. Also visible were π-cation electrostatic interactions with ARG-11 and LYS-183 at distances of 3.69 and 3.59 Å respectively, and π-anion electrostatic interactions with ASP-159 and ASP-196 at 3.76 and 3.73 Å respectively. The hydrophobic interactions include Amide π-stacked with ARG-195, alkyl interaction with LYS-199, and π-alkyl interactions with ARG-195, ILE-163, and CYS-192. The binding interactions of both compounds with 2HNL, unlike that with 1TU7, lacked electrostatic interactions.
Comparing the binding profiles 20 with FBZ revealed a significant similarity in their binding mode with both receptors.
Although both compounds showed more remarkable similarity in their binding pattern, their binding scores differ significantly. This may be attributed to the variations in the strength of their interactions and the contributions of Van der Waals interactions which were not visible.     Table VII. In addition, the 3D structural orientations of compound 20 in the active sites of both receptors before and after MD simulation were presented in Figures 18 and 19 for 1TU7 and 2HNL, respectively. Furthermore, the 2D views of the resulting binding interactions of the complexes before and after the MD simulation was shown in Figures 20 and 21 for 20_1TU7 and 20_2HNL, respectively. The overall average RMSD values were estimated as 9.4609 and 11.8643 Å for 20_1TU7 and 20_2HNL, respectively. Figure  15 shows that 20_2HNL deviated significantly (21.3988 Å) at the start of the simulation and fell gently as the simulation progressed, reaching equilibrium at 10 ns (0.0000 Å). For 20_1TU7, the RMSD increased gradually from 0.0000 Å at the start of the simulation and reached the peak (17.8253 Å) at 5.5 ns, after which it fell to 4.6270 Å at 10 ns. The RMSF is more like a calculation of the flexibility or the extent of movement of individual residue during a simulation. The RMSF plots in Figure  16 showed that the protein residues were slightly flexible between 0.5 and 1 Å during the trajectory, indicating the stable proteins during the simulation 30 . The Rg measures the degree of a protein's compactness during the trajectory. Decreasing Rg indicates reducing residues' flexibilities and more stability for the protein. The variations of Rg from the start to the end of the simulation were 0.60 and 0.52 Å for 1TU7 and 2HNL, respectively (Figure 17), connoting slight changes in the protein compactness as the simulation progresses, which therefore means both complexes could be stable. Furthermore, the result of binding free energy (MM/GBSA) estimated for 20_1TU7 and 20_2HNL (Table VII) shows that 20_1TU7 has a positive free energy value (+320.7898 kcal/mol) while 20_2HNL has a negative value of -29.6742 kcal/mol. This shows that the binding interactions of compound 20 with 2HNL are favorable, energetically stable, and bind strongly with the receptor. A similar observation was reported elsewhere for binding free energy change (MM/GBSA) calculated for some analogs of 2-aryl benzimidazole in a complex with PdxK 31 . Additionally, the orientations of compound 20 in the active sites of both receptors before and after the MD simulation was shown in Figures 18 to 19, which showed a slight change in the position of compound 20 on 2HNL after the simulation ( Figure 19B). However, the position of the compound on 1TU7 appeared to be reversed, curving away from the receptor's active site after the simulation (Figure 18B). Figures 20  and 21 presented the resulting binding interactions for 20_1TU7 and 20_2HNL, respectively. From Figure 20B, no interactions were visible between compound 20 and 1TU7, which could be attributed to the reversed position of the compound on the receptor after the MD simulation. From Figure 21B, on the other hand, none of the previous interactions was retained, though new interactions were formed, including one conventional Hydrogen bond with ARG-60 at a distance of 1.93 Å, and two carbon-hydrogen bonds with PHE-123 at 2.94 and 3.00 Å. Others are two electrostatic bonds, one π-cation with ARG-60 at 3.66 Å and one π-anion with PHE-225 at 4.40 Å, and one hydrophobic interaction (π-π stacked) with PHE-225 at a distance of 4.44 Å. Therefore, Prostaglandin D synthase (2HNL) is confirmed as a promising therapeutic target for benzoxaborole-benzimidazole hybrids.

CONCLUSION
The combined molecular docking screening and pharmacokinetics approach identified four benzimidazole-benzoxaborole analogs (14, 16, 19, and 20) as the most active molecules forming high binding affinities with the two protein receptors (2HNL and 1TU7). These molecules were predicted to be orally bioavailable, less toxic, and possessed better pharmacokinetic profiles than the reference compound (FBZ). These compounds' predicted pharmacological interaction profiles and FBZ generally fit well into the target site cavities. Compound 20 showed a typical characteristic of a multi-target drug molecule as it binds energetically stronger with both targets than other analogs. Additionally, the MD simulation revealed the stability of the interactions between compound 20 and 2HNL. Hence, the selected molecules, especially 20, could be developed and further evaluated as potential drug candidates for treating onchocerciasis. More so, FBZ remains promising and could still be improved to address the problem of oral bioavailability and the threats posed by its toxicity level.