Enrichments of Ensemble Docking Strategy Based on the Bayesian Model

Challenges remained in structure-based drug discovery which include protein exibility in binding site. Thus, concerning the exibility of proteins, docking into an ensemble of rigid conformations (ensemble docking) have been proposed with incorporation into protein exibility with expects that it could provide higher enrichments than rigid single receptor. Here we have developed the ensemble docking strategy by using Bayesian Model algorithms, and this method is validated by three proteins: BTK, JAK and PARP. The Bayesian Model was used to integrate independent docking runs of an ensemble of rigid crystal structures and MD simulations. Results The structure of MD simulations outperforms the crystal structures in separating inhibitors from decoys in BTK and PARP. Further, the results demonstrated that the ensemble docking strategy has better performance than rigid single conformation.


Introduction
As a computational methodology, structure-based drug design (SBDD) was introduced in drug discovery in the 1970s [1]. This predictive technique is valuable in limiting project timeline and reduces the cost of drug discovery. A large quantity of protein structures of drug targets was acquired with the blossom of experimental knowledge in crystallography, and thus structure-based drug design were applied in many cases to discover potential active candidates including antivirals for HIV and in uenza [2,3]. However, challenges remain in the area, such as the protein exibility in binding site, solvent treatment, and ligand electronic effects. Techniques that account for some degree of protein exibility are expected to increase the docking accuracy in both binding mode and a nity prediction [4].
Protein exibility could be accounted for through the simulation of receptor-ligand complex (induced-t) or an ensemble of rigid receptor conformations [5,6]. The rst time that ensemble docking of small molecules was taken into account as a strategy to accommodate the receptor exibility is conducted by multiple protein structures through exploring a number of conformational subsets of complex [7]. This subset of complex could be gathered from several ways including experimental data (X-ray crystallography and NMR spectroscopy) [8][9][10] and computational sampling (molecular dynamics simulations and homology models) [11,12].
Ensemble docking involves docking compounds into multiple conformations of a target rather than one single protein structure [13]. Some ensemble docking methods employ multiple receptor conformations in a single docking run during the pose-generation phase [14]. Others perform a series of independent docking runs based on all members of the ensemble followed by cluster analysis to select representative conformations [15,16]. The ensemble docking could yield higher enrichment of known inhibitor compared with single receptor conformation [17]. However, too many conformations used are unable to increase the accuracy. It is found that a 3-5 structures ensemble would be enough to yield satisfactory results more likely [9].
The key point is how to integrate the predictions of a series of independent docking runs based on multiple conformations for a speci ed target. A plausible way is merging the independent predictions into one entry and reranking or re-scoring. However, as the consensus scoring did, the strategy simply re-ranking and re-scoring using different docking algorithms is quite sensitive to each independent docking study [18]. An alternative way is the machine learning approaches that had been applied in virtual screening [19], which could integrate the predictions of different scoring functions for multiple conformations [20]. The Bayesian Model served as a support vector machine of machine learning algorithms has been a useful part of computer-aided drug discovery for many years and were popularized in Pipeline Pilot [21,22]. In 2004, it is observed that the enrichment factor of the consensus scoring had been improved by using Bayesian Model to integrate the predictions from ve scoring functions [23]. Besides, multiple receptors of ROCK1 had been combined as a single integrated model by using Bayesian Model, and the researchers thus found that the integrated model achieves the best prediction power [20,24].
In this study, we employed the Bayesian Model for integrating the predictions of molecular docking of each ensemble, evaluating and validating the performance of each single protein structure and ensemble of protein structures. The procedures of this ensemble docking via multiple receptor conformations are summarized in Fig. 1.

Methods
Considering the pharmaceutical importance of different targets in drug design/discovery, we have selected the Bruton's tyrosine kinase (BTK), Janus protein tyrosine kinases (JAK) and poly (ADP-ribose) polymerases (PARP) to develop an ensemble docking strategy. They are closely related to chronic immune diseases (rheumatoid arthritis and lupus) [25,26], cytokines intracellular signaling pathway [27], and cancer [28][29][30], respectively. In this work, we have employed GOLD docking program to perform molecular docking based on each crystal structure and MD (molecular dynamics) snapshot. Finally, an ensemble docking strategy was built to integrate the independent docking runs of each ensemble based on the Bayesian Model.

Ligand Dataset
The active ligands with quantitative biological activities (Ki, Kd, and IC50) were extracted from the Binding database [31]. The number of inhibitors of each target protein is 125, 152, and 134 for BTK, JAK, and PARP, respectively. The validation dataset was prepared by mixing the known inhibitors with twenty times of diverse decoys which were randomly selected from the ZINC database [32]. Thus, the number of molecules for BTK, JAK and PARP in validation dataset is 2625, 3192, and 2814, respectively.

Protein Structures Preparation
The corresponding crystal ligand-bound structures of target proteins were retrieved from RCSB Protein Data Bank (PDB, https://www.rcsb.org/). There 40, 22, and 38 crystal-structure complexes for BTK, JAK, and PARP, respectively.
The water and co-crystallized ligands of crystal structures were removed. Open-Source software PyMOL version 1.3 was employed to align structures of BTK, JAK, and PARP based on 4Z3V, 4EHZ and 6BHV, respectively.

Representative MD Simulations Generation
For each protein, the structure performing the best discrimination power was chosen as the initial structure of MD simulation. Representative MD simulations were determined from a trajectory using gromos algorithms [33] in GROMACS. The center of a cluster, the structure with the smallest average RMSD (root mean square deviation) from all other structures of the cluster, is determined as representative MD simulations. In detail, 2000 snapshots of the trajectory (20-40 ns) were extracted from the stable phase and RMSD matrix between every two protein conformations was built. The 2000 complexes were clustered into four categories when the RMSD cut off was set as 0.2 nm. For each category, the center of the cluster was selected as the representative structure. To evaluate the prediction power of each Bayesian Model, 70% of the inhibitors and non-inhibitors in the validation dataset were randomly selected for the training set, and the remaining 30% compounds for the test set (see Additional le 1).  Table S1), 6BBU (2.08 Å), 6N7A (1.33 Å), 6AAH (1.83 Å) and 4EHZ (2.17 Å) for JAK (Additional le 2: Table S2) and 4ZZZ (1.9 Å), 5WS1 (1.9 Å), 6I8T (2.1 Å) and 6NRH (1.5 Å) for PARP (Additional le 2: Table S3). The scatterplots of RMSD values for crystal structures versus residues were also calculated (Additional le 2: Figure S1).

Performance in Docking Calculations
We have previously studied that the compatible scoring functions by considering both reproducibility of ligand conformation and binding a nity prediction [34,35]. The scoring function CHEMPLP was demonstrated to be appropriate in both docking power and scoring power and thus applied to scoring function (Additional le 2: Table  S4).
The co-crystal ligand and prepared validation dataset was used to assess the docking power and discrimination power which presented by RMSD and P-value. The number of active compounds and decoys was listed in Table 1.
Docking calculations were conducted by screening the validation dataset to discriminate the inhibitors from noninhibitors. U-test (P-value) was used to evaluate the signi cance of the difference between the distributions of docking scores of the inhibitors and non-inhibitors.
The RMSD value below 2 Å is considered as a successful prediction of binding pose. Though, as shown in Table 2, in PARP the structure 6I8T is the unique that had succeed in pose prediction while four representative crystal structures of BTK and JAK have successfully reproduced the ligand conformation by re-docking the ligand as illustrated in Table 2 and Fig. 3. For BTK, each representative crystal structure also has The P-values suggested the signi cant difference of docking scores between inhibitors and decoys which means favorable discrimination power ( Table 2, Fig. 4). The best P-values of three targets were 8.586E-74, 2.438E-67 and 6.260E-57 for BTK, JAK, and PARP, respectively. Apparently, the crystal structures of BTK performed best in distinguish the inhibitors and non-inhibitors. The docking performance of different target were different, even the different structures for the same target exhibited remarkable difference. On the basis of P-values, the representative complex 5P9K, 6N7A and 6I8T was chosen as the initial structure of BTK, JAK and PARP for MD simulations, respectively.

Molecular Dynamics Simulations
The crystal structures 5P9K, 6N7A and 6I8T were chosen as the initial structure to construct the MD simulations. For each complex, 40 ns MD simulations were enough for trajectory to achieve equilibrium. The structural deviation between the structures in the dynamic simulation process and the initial structure were calculated based on the backbone. As shown in Fig. 5, the RMSD of backbone for BTK, JAK and PARP is under 0.3 nm which means the system has equilibrated in the molecule dynamic simulation process. To tolerate structural differences in the stable phase of each protein, 20-40 ns trajectory was collected. The RMSD matrix between each two conformations was constructed, based on which can perform the subsequent cluster analysis. On this basis, the 20-40 ns trajectory was clustered, and the conformation with middle RMSD obtained in the cluster can be used as a representative structure (Fig. 6). For BTK, the structures 5P9K_3768, 5P9K_2526, 5P9K_2992 and 5P9K_2629 were recovered as representative MD snapshots. While the structures 6N7A_3398, 6N7A_2014, 6N7A_3726 and 6N7A_2010 as well as 6I8T_3041, 6I8T_3906, 6I8T_2072 and 6I8T_3274 were recovered as representative MD simulations for JAK and PARP respectively.

The Correlation Analysis of Docking Performance Based on Crystal Structures and Structures from MD Simulations
Docking performance was evaluated by docking scores which were obtained from molecular docking simulation. The docking scores of the known inhibitors and non-inhibitors of each target based on each representative crystal structure and structure from MD simulations were used to execute the correlation analysis including Pearson and Spearman ranking correlation analysis. When absolute value of the correlation coe cient (r or ρ) is greater than or equal to 0.8, it is considered that the two variables are highly correlated. While when less than 0.3 indicates weak correlation, basically irrelevant. The redundant structures whose correlation coe cient is high ( 0.9), which means the compared two structures possess highly similar docking results, were discarded. Then the Pearson correlation coe cients (r) and the Spearman ranking correlation coe cients (ρ) of the docking scores based on each two compared protein structures were calculated (Fig. 7). In detail, the results were concluded in Additional le 2: Table  S5 and Additional le 2: The ROC curves were used to illustrate the virtual screening performance of Bayesian Models (Fig. 8), the corresponding AUC were displayed in Table 4. It is consistent with the Bayesian score; the MD simulations outperformed the crystal structures in BTK and PARP. However, the crystal structures combined ensemble (Ensemble 7) was presenting more outstanding performance than the MD simulations combined ensemble (Ensemble 8).
Lower correlation of docking scores between crystal structures compared with MD simulations may account for this. When different structures used for the docking calculations, the top ranking compounds are quite different. This elucidates the importance of choosing the appropriate protein structure in docking calculations, and then in ensemble constructing.

Conclusions
In this study, the ensemble docking was applied in BTK, JAK and PARP to retrieve known inhibitors from established ligand training set and test set and the performance were compared with standard single rigid receptor. By using Bayesian Model, the independent predictions of molecular docking of each ensemble were integrated. Additionally, the overall integrated predictions of molecular docking were represented by ROC curves to evaluate the enrichment of each single rigid receptor and ensemble. By applying this ensemble docking strategy in three targets, a total of 14 ensembles for BTK, JAK and PARP, were screened against the ligand dataset. The rigid crystal structures and MD snapshots were used in the same way to conduct the comparison between them. It is found that most of the ensembles are performing better than single rigid receptor in enrichment. Though, the ensembles in different size for three proteins suggest no absolute trend in enrichment of known inhibitors, in general, the enrichment is increasing with the ensemble size until 4-5 membered ensemble which also demonstrated by the academia before [9]. The larger size ensemble, such as six, seven or eight-size ensemble didn't express dominant superiority in enrichment. It suggests that the ensemble docking strategy built in this work could increase the enrichment of known inhibitors e ciently in BTK, JAK and PARP. We expect that the ensemble docking strategy can serve as a useful tool in virtual screening to nd more promising and diverse inhibitors when employing in more targets.

Declarations
Availability of data and materials The datasets supporting the conclusions of this article are included within the article and its additional le.

Competing interests
The authors declare that they have no competing interests.    The frequency distributions of GOLD docking scoring functions of actives and decoys for each of four representative protein structures which are crystallographically derived from of BTK, JAK and PARP, respectively.
Page 17/21 Figure 5 Time evolution of the RMSD of the receptor backbone compared with the initial structure for the representative complex 5P9K, 6N7A and 6I8T for BTK, JAK and PARP, respectively.  The heat maps were used to represent the signi cance of (left) Pearson correlation coe cients (r) and (right) Spearman ranking correlation coe cients (ρ) of the GOLD docking scores of the known inhibitors and noninhibitors for each pair structures of three targets.