{'fraction_non_alphanumeric': 0.03642243664510656, 'fraction_numerical': 0.03318842116424557, 'mean_word_length': 5.092700258397933, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
MOLECULAR MODELING AND DOCKING ANALYSIS OF BIS- INDOLYMETHANES DERIVATIVES AS HUMAN β-GLUCURONIDASE ENZYME INHIBITORS

Ibrahim* M. T. and Muhammad, U. Department of Chemistry, Faculty of Physical Science, Ahmadu Bello University, P.M.B 1045, Zaria, Kaduna State Nigeria. Department of Science Laboratory Technology, School of Technology, Kano State Polytechnic, Kano State Nigeria. *Corresponding author’s-mail: muhdtk1988@gmail.com, Phone Number: 08069651985 ABSTRACT β-glucuronidase enzyme is present mostly in mammals’ tissues. β-glucuronidase is present in kidney, bile, serum, urine and spleen. In eukaryotic and prokaryotic organisms, it is important in the process of breaking down of β-glucuronide. It also helps in the neutralization of reactivity of some metabolites that are associated to many diseases. The most stable geometry of the dataset were obtained adopting DFT method at B3LYP/6-31G* level of theory. The model was developed using MLR analysis adopting GFA method. Molecular docking was also performed to portray the binding mode of these bis-indolymethanes derivatives in the binding pocket of their target receptor (human βglucuronidase). The selected model was assessed and chosen based on its statistical fitness with Rtrng=0.907233, Radj=0.881465, Qcv=0.833795, and Rtest=0.609841.And also, the significance and impart of each physicochemical parameters to the selected model were determine by their ME values. Molecular docking analysis revealed that amino acid such asALA49, SER52, ASP53, PHE51, VAL96, LEU92, TYR188, TYR199 and PHE200 might be responsible for the most promised binding affinity of the reported docked ligands. The molecular docking results showed that the reported compounds were better than the standard β-glucuronidase inhibitor. The results of this findings paved way for designing novel β-glucuronidase inhibitors.


INTRODUCTION
β-glucuronidase enzyme is present mostly in mammals tissues, kidney, bile, serum, urine and spleen (Ali et al., 2016;Gloux et al., 2011). In eukaryotic and prokaryotic organisms, the enzyme is important in the process of breaking down of β-glucuronide (Beaud et al., 2005). It also helps in the neutralization of reactivity of some metabolites that are associated to many diseases (De Moreno de LeBlanc and Perdigón, 2005).It has been shown that increase in performance of this enzyme can lead to numerous unhealthy situations (Salar et al., 2016;Taha et al., 2015).This enzyme was stated to be sent to synovial fluid during inflammatory joint disorders (Taha et al., 2018).It is very paramount to devise a means to prevent the adverse effect of β-glucuronidase so as to stop many unhealthy situations caused by the enzyme. Due to their extensive uses in medicinal chemistry, pharmacology and biochemistry, bisindolymethanes were identified to possess different biological activities such antibacterial, HIV-1 integrase inhibitors, antitumor and antifungal, antimicrobial and aromatase inhibitors for breast cancer (Kamal et al., 2009;Lézé et al., 2004;Nagase et al., 2010).Also, some of these compounds are used by animals(humans) in the metabolism of estrogen to treat some sickness such as extended weakness, critical bowel symptom and fibromyalgia (Chakrabarty et al., 2002). Computational chemistry is a unique area in the drug design and development arena which provides in-silico methods and software that are employed in the discovery and production of new compounds of medicinal benefit (Jorgensen, 2004).
Quantitative structure-activity relationships (QSAR) is an in-silico method used to correlate the response variable (biological activities) with different descriptors (physicochemical properties) associated with the structures of a particular molecule (Ojha Lokendra et al., 2013).While an in-silico method used to predict the binding energy of http://dx.doi.org/10.4314/bajopas.v14i1.4 intermolecular complexes based on their 3D structures is known as molecular docking (Kitchen et al., 2004).This study is aimed at carrying out QSAR and molecular docking analysis on bis-indolymethanes derivatives against β-glucuronidase enzyme.

MATERIALS AND METHODS QSAR modelling methodology
Thirty two (32) derivatives of bisindolymethanes and their β-glucuronidase inhibitory activities (IC50) were retrieved from the work of Taha et al., (2018) for the purpose of this study. After data retrieval from their source, the inhibitory activities IC50 in (μM)of the studied data were transformed to their corresponding negative logarithm scale (pIC50)using equation 1 in order to increase linearity in the activities value. Chemdraw software was adopted for drawing the structures of all the studied data (Ibrahim et al., 2019). pIC50 = log (1/IC50)(1) In determining the structures of all the data set at global minima on Potential energy surface (PES) (stable structure), Density functional method (B3LYP/6-31G* level of theory) was employed to achieve the searching of the stable structures of all the dataset on potential energy surface (Amin and Gayen, 2016).For the generation of the physicochemical descriptors, the already optimized structures were save in SDF a file format recognized only by the Pharmaceutical data exploration laboratory tool kit (PaDEL descriptor tool kit).PaDEL descriptor tool kit was used to compute both 1D, 2D and 3D descriptors(Yap, 2011). Before data set splitting, the data were pre-treated using data pre-treatment software retrieved from drug theoretical and cheminformatics Laboratory(DTC Lab) to remove redundant and constant values from the data (Ambure et al., 2015). Data division software was further used to split the data into model building set (75%) and validation set (25%) (Kennard and Stone, 1969). The model building set was used to generate the models using multi-linear regression analysis adopting genetic function algorithm method. The equation for the regression analysis is shown in equation (2). Y = A1x1+ A2x2 + A3x3 + C(2) where Y is the pIC50(dependent variable), 'A's are coefficients for the descriptors(which are the 'x's), and 'C' is the constant for the regression equation (Ibrahim et al., 2020b). After generating the models, it is very important to assess the high predict power, reliability, stability and robustness of the generated models using the squared of the correlation coefficient (R 2 ), cross-validation coefficient (Qcv 2 ), and adjusted squared of the correlation coefficient (Radj 2 ) of the model (Jalali-Heravi and Kyani, 2004;Tropsha and Bajorath, 2015)the equations for these listed validation parameters are defined as: Qcv 2 = 1- where Ypred is the predicted pIC50,Yexp is the observed pIC50respectively of the validation set and Ymntrngis the averagepIC50 value of the model building set. Variation inflation factors (VIF)is also important in QSAR which is used to determine the multicollinearity problem of the physicochemical parameters(descriptors) in aQSAR model, If VIF values isone (1), there is no multicollinearity problem/inter-correlation between the variable. But if VIF values is between one (1) to five (5), the selected model can be accepted and therefore regard as valid and if VIF values is greater than ten (10), therefore the selected model is bad and therefore rejected (not free from multicollinearity problem/inter-correlation) (Beheshti et al., 2016). VIF can be determine using equation 6 below: VIF = 1 1− 2 (6) where R 2 is the correlation coefficient of the model. The mean effect (ME)is employed to determine the degree of contribution and significance of individual physicochemical descriptors to the selected model which indicates the direction in the activities of the compounds whether increase or decrease against their target enzyme. Mean effect help in ligand-based drug design by giving a hint on which physicochemical descriptor to give much consideration when carrying out structural modifications on the template. It is given by the expression below: where βj is the coefficient of the physicochemical parameter J in that selected model, dij is the value of the physicochemical parameter in the data matrix for each molecule in the model building set and MFj is the mean effect of physicochemical parameter j in the selected model, m is the number of physicochemical parameter that appear in the selected model and n is the number of molecules in the model building set (Ibrahim et al., 2020a). Domain of applicability is very important in QSAR model validation most especially in the quality of the model predictions and control of potential misuse of the models outcome. Also, it helps to figure out influential and outliers among the compounds in the data. The domain of applicability of the model must be exploited (Roy et al., 2017). As such leverage approach was adopted in this case and is given as hi: where Y is p × q independent variable matrix of the model building set compounds, yi is the model building compounds matrix I, and Y T is the transpose matrix Y utilized in developing the model. The thresh-hold value (h*)as indication tool and is the boundary for Y values and is given as: where z is the number of compounds in the model building set and qis the number of independent variable in the selected model. For any QSAR model to be considered as valid and used, it has to pass the Internal and external validations assessment (Veerasamy et al., 2011).

Molecular docking simulation methodology
Docking simulation was performed to study the nature and mode of binding interactions between the binding pocket of human βglucuronidase and the ligands utilizing Discovery studio visualizer, Autodock Vina of Pyrex virtual screening and UCSF Chimera docking software. The coordinates and dimensions of the grid box used for the docking simulation are X: 81.5147 Å, Y:90.5618Å and Z:138.5886Å respectively. Ligands were prepared prior to the commencement of the docking simulation, by saving the optimum conformation ascertained using density functional theory in protein data bank file (pdb file format). The crystal structure of Humanβ-glucuronidase was retrieved from pdb with pdbID 1bhg (Ibrahim et al., 2020c). The preparation of the human β-glucuronidase for the docking simulation was done using Discovery Studio Visualizer, by removing chain B, heteroatoms and co-ligands from the dimer saved also as protein data bank file (pdb file format) (Abdulfatai et al., 2017).Pyrex software was used in the execution of the docking simulation in which the ligands were docked to the binding site of the human β-glucuronidase (Trott and Olson, 2010). The complexes were rebuild using UCSF Chimera software for further investigation. The nature and mode of binding interactions of the complexes was investigated using the Discovery studio visualize (Abdulfatai et al., 2019).

RESULTS AND DISCUSSION QSAR modelling results
Four QSAR models were developed out of which the best model was selected and reported based its statistical significance. Model 1 was selected and reported as the best because of its statistical fitness. On comparing the statistical parameters of the selected model with those reported by Veerasamy et al., (2011) it can be seen that the statistical parameters of the selected and reported model were all greater than the minimum recommended values which confirmed the reliability of the model (Veerasamy et al., 2011). The squared correlation coefficient (R 2 trng) of the reported models was 0.907233 which means that the model can be able to explain about 90.72 % of the variations in the activities of these β-glucuronidase inhibitors (Golbraikh and Tropsha, 2002). Also the value of this R 2 trng(0.907233) was greater than that of its corresponding R 2 adjvalue (0.881465) which confirm the significance of the reported model (Ambure et al., 2015). The reliability of the reported model was further confirmed by the calculation of the predicted activities of the validation set compounds and the external validationR 2 testvalue (0.609841) (Roy et al., 2016 To confirm the quality of the selected model, the Predicted activities of both the model building set and that of the validation set were plotted against the actual activities (Figure 1). The indicator used in this case is R 2 value of both the plot and that of the internal validation, the quality of the selected model was confirmed by the corroboration of R 2 value (0.9072) of the plot and that of the internal validation (R 2 trng=0.907233). Every Good QSAR model is expected to be free from methodological/systematic error. In order to determine whether the selected and reported model is free from systematic error, the predicted activities were plotted against their standardized residuals ( Figure 2). The selected model was confirmed to be free from methodological error by even distribution of the standardized residuals on the plot. The observed activity was seen to have good correlations with the predicted activity. Table 1   The correlation analysis on the independent variables inthe model building set of the selected model in Table 2 indicates the importance of the independent variables to the model. The independent variables were found to have no correlation with one another as no two descriptors have their values close to one. The computed Variation Inflation Factor values for all the independent variables were found to be less than 5 (see Table 2) indicating the statistical fitness of the selected model and no multicollinearity problem exist between the independent variables. The mean effect(ME) value (Table 2) shows the degree of contribution of an independent variable, in comparison to others in the reported model. The positive or negative coefficients of the independent variable show the direction of the activity in inhibiting the β-glucuronidase enzyme whether high or low. From the mean effect, GATS2e(Geary autocorrelation -lag 2 / weighted by Sanderson electronegativities) gave the minimum degree of contribution with the negative value of -0.08929 which indicates that this physicochemical parameter contributes negatively to the potency of bis-indolymethanes against their target enzyme (β-glucuronidase) in the sense that if the number of this physicochemical parameter is reduced, it means that the potency of bis-indolymethanes will be high against β-glucuronidase and vice versa. On the other hand, the mean effect values forGATS3e(Geary autocorrelation -lag 3 / weighted by Sanderson electronegativities), GATS4s(Geary autocorrelation -lag 4 / weighted by I-state), SpMAD_Dzs(Spectral mean absolute deviation from Barysz matrix / weighted by I-state ) and SpMax5_Bhs(Largest absolute eigenvalue of Burden modified matrixn 5 / weighted by relative I-state)signifies their positive contributions toward the effectiveness of bis-indolymethanes against β-glucuronidase each with positive value of +0.219708, +0.105168, +0.200507 and +0.563911 respectively. It indicates that addition of these descriptors to the bis-indolymethanes will increase their potency against β-glucuronidase and vice versa. The trend in the individual contribution given by these descriptors is given as SpMax5_Bhs ˃ GATS3e ˃ SpMAD_Dzs ˃ GATS4s ˃ GATS2e The plot of leverage values calculated for all the dataset and the standardized residuals (Williams plot) (see figure 3), which permit a graphical identification of both influential and outliers compounds in the selected model (Beheshti et al., 2016). From the plot, all the compounds of the model building set and 4 from the validation set were within the domain of the model. And only four influential compounds were observed from the validation set. Those influential compounds can be said to have their mechanism of action different from those within the domain of applicability of the reported model. More so, there were no outliers in both model building set and the validation set with their standardized residual greater than the +3 or -3 standard deviation unit.

Results of Molecular docking analysis
Molecular docking simulation on all the thirty two (32) bis-indolymethanes derivatives was performed to investigate the mode of binding interactions between them and their target enzyme (human β-glucuronidase, pdb ID:1bhg). The binding energy of all the studied ligands ranges from -10 kcal/mol to -12.4 kcal/mol as shown in Table 1.

Trng Set Test Set
The next to ligand 28 reported with higher binding affinity is ligand 19 with -12.1kcal/mol binding affinity, it interacted with the human βglucuronidase through hydrophobic interactions  with amino acid chains ALA49, SER52, ASP53,  PHE51, SER52, VAL201, VAL96, LEU92, TYR188,  TYR199 and PHE200 and also via conventional  hydrogen  bond  interactions with:HIS94 (2.10014) and PHE51 (2.53095) amino acid residues. Also, Ligand 1also shows good interaction with high binding affinity of -11.8kcal/mol. Hydrophobic interactions with GLN202,TRP90,ALA49,PHE51 and PHE95, Electrostatic interactions with ALA49, SER52, ASP53, PHE51, SER52, VAL96, LEU92, TYR188, TYR199 and PHE200, and conventional hydrogen bond interactions withPHE200 (2.97636) and PHE51 (2.63555) were observed with mentioned amino acid residues of the human β-glucuronidase. Beside the mentioned ligands, ligand 25 with binding affinity of -11.7 kcal/mol was also observed to interact with the binding pocket of the human β-glucuronidase through hydrogen, hydrophobic and electrostatic interaction as shown in Table 3.The standard drug (D-saccharic acid 1,4-lactone) with the binding affinity of -5.7kcal/mol formed only hydrogen bond with ASN502(2.42729Å) and GLN524 (3.44886Å)amino acid residues of the human β-glucuronidase. All the compounds were seen to be more active than the standard drug. Figure 4 shows the 3D structure of the reported compounds investigated using PyMOL.  CONCLUSION QSAR modelling on some bis-indolymethanes was conducted using Genetic Function Algorithm (GFA). The most stable geometry of the studied data were obtained using DFT method utilizing B3LYP/6-31G* level of theory.The selected model was assessed and chosen based on its statistical fitness with R 2 trng=0.907233, R 2 adj=0.881465, Qcv 2 =0.833795, and R 2 test=0.609841.Molecular Docking simulation reported between some selected compounds (compound 28, 19,1 and 25) and binding site of human β-glucuronidase enzyme, showed that these amino acid residuesALA49, SER52, ASP53, PHE51, VAL96, LEU92, TYR188, TYR199 and PHE200might be responsible for the most promised binding energy ofthe reported docked ligands.The reported compounds werefound to be more active than the standard drug used as control in this study. The result of this in-silico findings paved way for designing new novel βglucuronidase inhibitors. Author's contribution: Ibrahim M. T. and Muhammad Umma: Conducted the research Conflict of interest: Authors declare no conflict of interest.