Fragmented blind docking: a novel protein–ligand binding prediction protocol

Abstract

In the present paper we propose a novel blind docking protocol based on Autodock-Vina. The developed docking protocol can provide binding site identification and binding pose prediction at the same time, by a systematical exploration of the protein volume performed with several preliminary docking calculations. In our opinion, this protocol can be successfully applied during the first steps of the virtual screening pipeline, because it provides binding site identification and binding pose prediction at the same time without visual evaluation of the binding site. After the binding pose prediction, MM/GBSA re-scoring rescoring procedures has been applied to improve the accuracy of the protein–ligand bound state. The FRAD protocol has been tested on 116 protein–ligand complexes of the Heat Shock Protein 90 – alpha, on 176 of Human Immunodeficiency virus protease 1, and on more than 100 protein–ligand system taken from the PDBbind dataset. Overall, the FRAD approach combined to MM/GBSA re-scoring can be considered as a powerful tool to increase the accuracy and efficiency with respect to other standard docking approaches when the ligand-binding site is unknown.

Communicated by Ramaswamy H. Sarma

Keywords:

• Blind docking
• molecular docking
• protein–ligand interactions
• virtual screening

1. Introduction

Discovery and development of a new drug is a complex process which requires a lot of time and resources. The main goal of drug discovery is to obtain compounds that successfully interact with therapeutic targets (Batool et al., 2019; Hughes et al., 2011). Computer Aided Drug Design (CADD) approaches are widely used to increase the efficiency of the drug design (Macalino et al., 2015). CADD techniques represents a fundamental tool in the preliminary stages of drug discovery pipeline (Surabhi & Singh, 2018; Veselovsky & Ivanov, 2003). Computational drug design tools can be divided into ligand-based drug design (LBDD) and structure-based drug design (SBDD). The choice between these two methods dependents on the available information about the biological target (Aparoy et al., 2012). SBDD approach uses 3 D structure of the protein target for the screening of potential ligands (Batool et al., 2019). On the other side, LBDD approach is employed where 3 D structure of the protein–ligand complex is missing, and computational modelling can be applied to develop theoretical predictive models starting from the chemical description of the ligands with known biological activity (Bacilieri & Moro, 2006; Śledź & Caflisch, 2018). In this framework, molecular docking is one of the most powerful techniques of SBDD (Gupta et al., 2018). The aim of molecular docking is to predict the experimental binding mode of a small molecule within the binding site of the selected target receptor (Ferreira et al., 2015; Hernndez-Santoyo et al., 2013; Preto et al., 2018). A search algorithm and an energy scoring function are the basic tools of a docking methodology for generating and evaluating the ligand conformations (Li et al., 2019; Novič et al., 2016; Prieto-Martínez et al., 2018). The docking procedure has already demonstrated its ability to correctly predict the ligand binding pose starting from the preliminary knowledge of the binding region (Hetényi & Van Der Spoel, 2006; Seeliger & De Groot, 2010). In this case, a small docking box is selected around the protein binding site in order to facilitate the docking by focusing sampling of the translational, rotational, and torsional degrees of freedom of the ligand (Morris et al., 2008). However, there are situation in which the information about binding site is missing and it becomes necessary to explore the entire protein surface by docking algorithms (Hetényi & Van Der Spoel, 2006). Several methods have been developed to address this issue (Ghersi & Sanchez, 2009; Grosdidier et al., 2009; Hassan et al., 2017). AutoDock Vina (Vina) (Trott & Olson, 2010) is a widely employed technique when information about the binding site is available, ensuring a successful docking procedure. When the protein binding pocket is unidentified, Vina can execute the so-called blind docking (BD). In BD, the target is included into a single simulation box covering the entire surface of the protein (Muscat et al., 2020). This method has many limitations because it is difficult to exhaustively sample the whole energy landscape in a fixed number of steps to find the best ligand conformation (Hassan et al., 2017). Another approach to address the above-mentioned issue consists in reducing the search space, focusing only in some areas of the protein. The SiteHound algorithm (Hernandez et al., 2009) predicts the location of potential binding site and after that, and multiple independent docking procedures on smaller boxes centered on predicted binding sites can be performed. The results show that the docking focused on a small number of predicted binding sites reduces the computational time required to obtain the solution and generates more accurate results in terms of correct pose prediction in comparison to BD. However, if the real binding site is not included in the boxes in which the docking is carried out, the procedure could lead to an incorrect result. The main limitation of the two analyzed methods consists in a too inaccurate sampling box for the BD method and uncertain prediction of the binding pocket for the SiteHound algorithm. In this work, a fragmented docking protocol (FRAD) has been developed to improve the performance of the blind docking performed by Autodock-Vina. The main idea is to divide the protein volume in multiple overlapping boxes. A docking calculation has been performed for each of the previously-mentioned boxes, in order to cover the entire protein surface. The predicted binding poses have been collected, obtaining an ensemble of several protein–ligand binding conformations. Then, the protein–ligand bound states have been ranked considering the corresponding energy scores estimated by Autodock-Vina. The partition in several boxes allows the systematical and exhaustive exploration of the whole protein surface, which improves the discovery of ligand conformations adopted within each examined box. In addition, the complete investigation of the whole protein structure, intrinsically leads to discover of the real binding site. The FRAD molecular docking has been applied to 116 crystal structures of Heat Shock Protein 90 – alpha (Hsp90α), 176 of Human Immunodeficiency virus protease 1 (HIV-1 PR), and more than 150 protein–ligand complex coming from the PDBbind dataset (Z. Liu et al., 2015; R. Wang et al., 2005 ).

2. Material and methods

2.1. Proteins dataset

Two proteins have been considered for testing the protocol: The Heat Shock Protein 90 – alpha (Hsp90α) and the Human Immunodeficiency virus protease 1 (HIV-1 PR). The main reason for the choice of the latter proteins is because they are among the most widely used proteins in literature for docking purposes (Ahmed et al., 2013; Al-Sha’er & Taha, 2012; Chang et al., 2010; Mey et al., 2016; Mobaraki et al., 2019; Penkler & Tastan Bishop, 2019).

Heat Shock Protein is a globular chaperone protein that assists other proteins to fold correctly and allows cells to survive in extremely heat conditions. It can preserve the integrity of the cells if they are exposed to high temperatures by regulating the flux of calcium ions, maintaining the chromosomal stability and safeguarding the endoplasmic reticulum proteins homeostasis. If on the one hand, Hsp90 aid the human body to survive at elevate temperature, on the other hand, can stabilize the proteins necessary for tumor growth. For this reason, Hsp90 inhibitors are mainly investigated as anti-cancer drugs (Prodromou et al., 1997). Hsp90α is an isoform of Hsp and it is located in the cytoplasm of the eukaryotic cells. It consists of four structural domains: N-terminal domain, middle region, C-terminal domain and a linker region that connect the first two domains. The binding pocket is situated in the N-terminal domain, which is predominantly constituted by hydrophobic residues. The pocket is a binding site for many molecules included antibiotics like radicicol and geldanamycin, which have anti-tumor activity and ATP (Garrido et al., 2001; Prodromou & Pearl, 2003).

HIV-1 PR is a retroviral aspartyl protease, an enzyme which acts in peptide bond hydrolysis in retroviruses, which is necessary for the life cycle of the retrovirus HIV that causes AIDS. HIV-1 PR is constituted from 99 amino acid and exists as a homodimer with only one active site. The binding pocket is located between the identical subunits. Each monomer consist of a wide β-sheet region (a loop of glycine) which partly constitutes the binding site of the substrate and one of the two fundamental residues of aspartyl, Asp-25 and Asp-25′ which are at the bottom of the cavity (Brik & Wong, 2003; Deeks et al., 1997). The HIV-1 PR enzyme activity can be inhibited from HIV protease inhibitors by blocking the active site of the protease.

The docking protocol has been also tested on 100 protein–ligand complex coming from the PDBbind dataset (Z. Liu et al., 2015; R. Wang et al., 2005).

2.2. Dataset Preparation

FD and BD methods have been carried out on the same set of complexes extracted from the Protein Data Bank (PDB) (Berman et al., 2000). The dataset of Hsp90α is composed of 116 crystal structures, while the dataset of HIV-1 PR is composed of 176 crystal structures. As normal procedure, all the waters molecules have been removed from each PDB entry. Missing residues and atoms have been reconstructed to the main protein chain utilizing the program Modeller (Webb & Sali, 2016). The particular cases where the protein or ligand have a PDB record of double spatial positions of the atoms have been treated selecting only one of the positions pair. Hydrogens and Gasteiger charges have been added both the ligand and the protein employing python scripts from MGLTools (Morris et al., 2009; Sanner, 1999).

2.3. Docking procedure

2.3.1. FRAD protocol

In the FRAD protocol, the protein volume has been divided in multiple overlapping boxes, as shown in Figure 1. The size of each box has been chosen as twice the ligand size and the overlap between different boxes has been set to fifty percent. A docking calculation has been performed by Autodock-Vina (Trott & Olson, 2010) for each of the previously-mentioned boxes, in order to cover the entire protein surface. The exhaustiveness value has been set to 8. The first ten predicted binding poses have been collected, obtaining an ensemble of several protein–ligand binding conformations. Then, the protein–ligand bound states have been ranked considering the corresponding energy scores estimated by Autodock-Vina.

Figure 1. Main idea of FD: a docking procedure is performed in each box which moves on the whole protein surface. Protein is represented in blue, ligand in yellow.

2.3.2. Blind docking protocol

A single docking calculation has been performed by using Autodock-Vina (Trott & Olson, 2010) on the whole protein volume. The exhaustiveness value has been set to 8. Results of FRAD and BD methods have been compared by root mean squared deviation (RMSD) of ligand heavy atoms of the solution for FRAD and BD with respect to the experimental ligand structures. In literature a value of RMSD lower or equal to 0.2 nm is recommended as a limit value for a good pose reproduction (Cole et al., 2005; Gabb et al., 1997; Lape et al., 2010; Nissink et al., 2002; Zaheer-ul-Haq et al., 2010). The RMSD comparison has been made employing the ligand conformation to which Vina attributed as the best in the affinity value ranking.

2.4. MM-GBSA rescoring

Docking algorithms generate binding poses of compounds and rank them by means of scoring functions (Du et al., 2016; Kairys et al., 2019; Vuignier et al., 2010). However, docking scores and experimental binding affinities usually do not correlate, because screening large numbers of compounds in a reasonable time requires the use of approximate scoring function. Hence, docking results can be improved employing post-processing strategies (Quiroga & Villarreal, 2016; Rastelli & Pinzi, 2019). The MM-GBSA methods has been successfully applied in literature as re-scoring procedure after molecular docking (Chen et al., 2016; Genheden & Ryde, 2015; Sun, Li, Shen, et al., 2014; Zhang et al., 2017). The MM-GBSA parameters were set accordingly to previous literature (Muscat et al., 2020; Su et al., 2015). In detail, the bondi, mbondi, and mbondi2 radii sets were prepared using the antechamber program of AmberTools19 (Salomon-Ferrer et al., 2013). The optimized version GBneck (GBneck2) (Nguyen et al., 2013) was used as GB model. The default setting of MM-GBSA surface tension (α = 0.0072 kcal/mol Å2) and the non-polar free energy correction term (β = 0) were applied. The dielectric constant of the medium was set to 80. Successful rescoring is obtained if the lowest binding free energy value belongs to the pose with a RMSD less than 0.2 nm.

3. Results

In this section, results regarding the FRAD and BD docking procedures are reported to compare performances of both methods as applied to ligand pose discovery in Hsp90α and in HIV-1 PR datasets. As a first result, we’ve shown if the binding poses have been correctly predicted. Then if the rescoring method produces improvements to the results. RMSD Tables concerning the use of FRAD and BD on 116 Hsp90α and 176 HIV-1 PR proteins with co-crystallized ligands are reported in Supplementary Information (Tables S1 and S2). The FRAD protocol has been also tested on 100 protein–ligand systems taken from the PDBbind dataset (Z. Liu et al., 2015; R. Wang et al., 2005), as reported in Table S3.

3.1. FRAD vs BD in protein–ligand complexes of Hsp90α

As described in the previous sections, several binding ligand poses have been generated for each protein and the pose considered for the validation is the one which Vina has attributed the maximum affinity value. The example shows in Figure 2(A) qualitatively shows the FRAD ability to produce more accurate poses with respect to the BD method. Moreover, Figure 2 exhibits a comparison between the protein–ligand interactions in crystallized complex (Figure 2(B)), in complex generated employing the FRAD method (Figure 2(C)) and in complex generated by the BD method (Figure 2(D)). It is worth mentioning that FRAD approach can identify the exact number of h-bonds with the correct amino acids and is able to conserve a high fidelity of the hydrophobic contacts (Figure 2(C)), if compared with BD method (Figure 2(D)).

Figure 2. Example of docking result employing FRAD and BD methods (PDB: 2yki). (A) Protein is represented in blue, the ligand in the crystal structures in red, the ligand conformation generated with FRAD method in cyan (RMSD is 0.04 nm) and the ligand conformation generated with BD method in orange (RMSD is 1.55 nm). (B) Protein–ligand interactions in crystallized complex. (C) Protein–ligand interactions in the generated with FRAD method complex. (D) Protein–ligand interactions in the generated with BD method complex. Elements in figure (A), (B) and (C) are the conformation of the ligands, hydrogen bonds in green and hydrophobic residue names.

In order to have a clear picture of the FRAD and BD performance, pie chart is selected as graphical tool to highlight the accuracy of both methods. The pie chart in Figure 3. Percentage of poses generated with RMSD lesser and greater than 0.2 nm compared the crystal structure pose using the FRAD and BD method. Results of FRAD methods are shown in the pie chart on the left and results of BD method are shown in the pie chart on the right. Red and blue indicate the poses with RMSD lesser and greater than 0.2 nm, respectively. (A) The comparison has been made considering the poses with the highest affinity value. An accurate ligand conformation has been discovered in 62.1% of cases using FRAD method and 27.6% using BD method. (B) The comparison has been made considering the poses with the best RMSD value between first 10 top poses. An accurate ligand conformation has been discovered in 91.4% of cases using FRAD method and 43.1% using BD method. A shows the percentages of generated ligand poses with RMSD value lesser and greater than 0.2 nm in comparison to the experimental ligand conformation. Considering the output conformation with greatest affinity, BD finds the correct pose prediction in only 26.7% of attempts, whereas FRAD finds the correct prediction in the 62.1% of attempts. The ligand conformation with the greatest affinity identified by the BD approach is not always the one closest to the real pose of the ligand. In this framework, the first ten conformations have been selected and the RMSD analysis between them and the original ligand has been computed. The hypothesis is to find the pose with the lowest RMSD value between the top ten records of Vina’s output conformations in order to evaluate if there is a pose predicted correctly with an RMSD less than 0.2 nm. Figure 3 shows the pie chart percentages of the obtained results considering the pose with lowest RMSD between the first ten poses. Interestingly, an improvement in the performance is notable for both FRAD and BD methods. The results highlight that the correct pose is generated among the ten top poses and the correct identification of it can be achieved through rescoring methodologies. In detail, FRAD finds in 91.4% of times the poses with RMSD value lower than 0.2 nm, while BD discovers in 43.1% of times the correct conformation.

Figure 3. Percentage of poses generated with RMSD lesser and greater than 0.2 nm compared the crystal structure pose using the FRAD and BD method. Results of FRAD methods are shown in the pie chart on the left and results of BD method are shown in the pie chart on the right. Red and blue indicate the poses with RMSD lesser and greater than 0.2 nm, respectively. (A) The comparison has been made considering the poses with the highest affinity value. An accurate ligand conformation has been discovered in 62.1% of cases using FRAD method and 27.6% using BD method. (B) The comparison has been made considering the poses with the best RMSD value between first 10 top poses. An accurate ligand conformation has been discovered in 91.4% of cases using FRAD method and 43.1% using BD method.

In Supplementary Information, Figure S1 shows the relationship between the improving in performance with the increase in the number of considered poses.

3.2. FRAD vs BD in protein–ligand complexes of HIV-1 PR

The second structure considered for the Vina FRAD and BD method is HIV-1 PR protein. In this context, the docking procedure has been performed 176 different ligand-protein complexes. As applied before for the Hsp90α protein, several binding ligand poses have been generated for each protein and the pose considered for the validation is the one which Vina has attributed the maximum affinity value with the protein. Figure 4(A) shows how an accurate pose has been predicted using FRAD method compared to the BD method. Moreover, Figure 4 exhibits a comparison between the protein–ligand interactions in crystallized complex (Figure 4(B)), in complex generated employing the FRAD method (Figure 4(C)) and in complex generated employing BD method (Figure 4(D)). It is worth mentioning that FRAD method is able to identify the exact number of h-bonds with the correct amino acids and is able to conserve a high fidelity of the hydrophobic contacts (Figure 4(C)), if compared with BD method (Figure 4(D)).

Figure 4. Example of docking result employing FRAD and BD methods (PDB: 1msn). (A) Protein is represented in blue, the ligand in the crystal structures in red, the ligand conformation generated with FRAD method in cyan (RMSD is 0.03 nm) and the ligand conformation generated with BD method in orange (RMSD is 2.08 nm). (B) Protein–ligand interactions in crystallized complex. (C) Protein–ligand interactions in the generated with FRAD method complex. (D) Protein–ligand interactions in the generated with BD method complex. Elements in figure (A), (B) and (C) are the conformation of the ligands, hydrogen bonds in green and hydrophobic residue names.

The pie chart in Figure 5 shows the percentages of generated ligand poses with RMSD value lesser and greater than 0.2 nm in comparison to the crystal structure ligand conformation. Considering the output conformations with greatest affinity, BD finds the correct pose prediction in only 31.3% of attempts, while FRAD finds the correct prediction in the 64.8% of attempts. Considering also this second dataset, the ligand conformation which Vina assigns the greatest affinity is not always the one closest to the real pose of the ligand. In this framework, the first ten conformations have been selected and the RMSD analysis between them and the experimental ligand has been computed. The hypothesis is to find the pose with the lowest RMSD value between the top ten records of Vina’s output conformations in order to evaluate if there is a pose predicted correctly with an RMSD less than 0.2 nm. Figure 5 shows the pie chart percentages of the obtained results considering the pose with lowest RMSD between the first ten poses. As pointed out in the previous section, an improvement in the performance is notable for both FRAD and BD methods. The results highlight that the correct pose is generated among the ten top poses and the correct identification of it can be achieved through rescoring methodologies. In detail, FRAD finds in 84.7% of times the poses with RMSD value lower than 0.2 nm, while BD discovers in 41.5%. of times the correct conformation.

Figure 5. Percentage of poses generated with RMSD lesser and greater than 0.2 nm compared the crystal structure pose using the FRAD and BD method. Results of FRAD methods are shown in the pie chart on the left and results of BD method are shown in the pie chart on the right. Red and blue indicate the poses with RMSD lesser and greater than 0.2 nm, respectively. (A) The comparison has been made considering the poses with the highest affinity value. An accurate ligand conformation has been discovered in 64.8% of cases using FRAD method and 31.3% using BD method. (B) The comparison has been made considering the poses with the best RMSD value between first 10 top poses. An accurate ligand conformation has been discovered in 84.7% of cases using FRAD method and 41.5% using BD method.

In Supplementary Information, Figure S2 exhibits the relationship between the improving in performance with the increase in the number of poses considered. A rescoring procedure has been performed in order to evaluate if MM-GBSA is able to attribute the best affinity value to the conformation with the lowest RMSD value, namely, to further validate and improve the FRAD method.

3.3. MM-GBSA rescoring

The evaluation of Vina's pose made considering the RMSD as principal parameter for the ranking is a commonly used protocol (Castro-Alvarez et al., 2017; Cavasotto & Abagyan, 2004; Hu et al., 2004; Prieto-Martínez et al., 2018), but the pose with the lowest RMSD does not always correspond to the one with greater affinity. The rescoring procedure may have a strong influence in the BD protocol, since the FRAD protocol is supposed to be performed within a variable protein region. To address this issue, a rescoring protocol on outgoing poses has been accomplished. Binding free energy has been calculated using the MM-GBSA method on the first ten protein–ligand complex whose affinity is the highest calculated by Vina. Ten protein–ligand complexes has been considered for the calculation. In Figure 6 the rescoring results for Hsp90α and HIV-1 PR structures are represented. In 75.7% of cases the greatest binding affinity estimated with MM-GBSA corresponds to ligand pose with a RMSD lesser than 0.2 nm compared the ligand in the crystal structure. The percentage of correct poses prediction is higher comparing the results in which the ranking is done with the affinity estimated by Vina. In HIV-1 PR the performance is about the same as that obtained by considering the docking affinities.

Figure 6. Rescoring results. The percentage of the poses generated with RMSD lesser and greater than 0.2 nm compared the crystal structure pose using FRAD method are shown for Hsp90α and HIV-1 PR structures. Results of Hsp90α are shown in the pie chart on the left and results of HIV-1 PR are shown in the pie chart on the right. Blue and orange indicate the poses with RMSD lesser and greater than 0.2 nm, respectively. An accurate ligand conformation has been discovered in 75.7% of cases in Hsp90α and 62.9% in HIV-1 PR.

4. Discussion

Protein–ligand docking is a powerful tool in drug discovery to predict binding modes and affinities of ligand (Ben-Shimon & Niv, 2015; Hetényi & van der Spoel, 2002; Hetényi & Van Der Spoel, 2006). The blind docking is a common strategy employed when the binding site of a target is unknown. However, BD has many limitations because it is unlikely to exhaustively sample the whole energy landscape to find the correct ligand pose (Hetényi & Van Der Spoel, 2011). One of the problems which may lead in a lack of accuracy in the blind docking results is due to an incorrect identification of the binding site (Brown & Jagt, 2004; Ghersi & Sanchez, 2009). Nowadays, the research community involves conspicuous effort in developing new software, like METADOCK2 (Imbernón et al., 2021), which increases the blind docking accuracy and efficiency in comparison with the Vina’s protocol. The speed up of this approach is due to the adoption of metaheuristic schemes. The strategies proposed usually employ specific software to find the active site in the target and then docking the ligands into the discovered binding site (Blaszczyk et al., 2016; Ghersi & Sanchez, 2009; Lin et al., 2016; Y. Liu et al., 2020). In particular, algorithms like SiteHound (Hernandez et al., 2009), CurPocket (Y. Liu et al., 2020) and COACH (Yang et al., 2013) have allowed an increased docking accuracy due to a more precise identification of the ligand binding site. These approaches have demonstrated how the blind docking results without information of the correct binding site are characterized by a very low success probability. However, if the true binding site is not included in the boxes in which the docking is carried out, the procedure will lead to an incorrect result. Another possibility is to perform a preliminary BD calculation to predict the binding site, and then a re-docking focused in the correct protein region (Agrawal et al., 2019). Despite a slight increase in accuracy, the computation cost of this procedure is higher and does not fully address the correct binding site prediction. The main goal of this work is to improve the BD procedure through a docking protocol called FRAD. In the FRAD protocol, the protein volume has been divided in multiple overlapping boxes. A docking calculation has been performed by Autodock-Vina (Trott & Olson, 2010) for each of the previously-mentioned boxes, in order to cover the entire protein surface. The first predicted binding poses have been collected, obtaining an ensemble of several protein–ligand binding conformations. Then, the protein–ligand bound states have been ranked considering the corresponding energy scores estimated by Autodock-Vina. A first comparison has been carried out between the RMSD values obtained considering only the ligand pose with the lowest docking energy. In this context, the performances obtained with the FRAD method are better than the BD method in both protein–ligand complexes analyzed. From another point of view, both the FRAD and BD methods increased their performance if we take into consideration the first ten poses, crowing nevertheless FRAD algorithm as the best one even in this case. The method of performing the docking procedure in several boxes that translate on the protein has proven to be a winning strategy in finding the lowest RMSD (<0.2 nm) conformation of the protein–ligand complex, namely the closest to the crystal structure. The improvement induced by this method has been confirmed by comparing the performances obtained in other works in which a BD has been carried on Hsp90α and HIV1-PR (Mey et al., 2016; Mobley et al., 2014).

A second analysis has been accomplished in order to evaluate the performance of Vina's scoring function, as done before in literature (Gaillard, 2018; Quiroga & Villarreal, 2016). Docking algorithms, like Vina, have been developed with the goal to screening large numbers of compounds in a reasonable time, for this reason, they use approximate scoring functions (Graves et al., 2008). It is reasonable to hypothyze that the best pose, namely the one with the minimum RMSD from the crystal structure, is not always the first in the rank. Indeed, it has become general opinion that docking results should be improved by means of more rigorous post-docking processing strategies (Palacio-Rodríguez et al., 2019; Quiroga & Villarreal, 2016; Rastelli & Pinzi, 2019). Several works have employed re-scoring approaches to predict binding affinities and afterwards compare these results with the docking scores (Lyne et al., 2006; Sgobba et al., 2012; Thompson et al., 2008). In this context, a study on Protein kinase proved that the correlation between calculated binding free energies with MM-GBSA and MM-PBSA methods and experimental values is higher than using docking scores (Sun, Li, Tian, et al., 2014). Instead, in a set of β-Amyloid Cleaving Enzyme 1 (BACE-1) the re-scoring docking poses using MM-GBSA did not improve the correlation with experimental affinities due to the ligand dataset (El Khoury et al., 2019). It is important to mention that BACE-1 ligands dataset consisted of macrocycles which have multiple flexible bonds that generate a large conformational space and for this reason require a more accurate MM-GBSA protocol(El Khoury et al., 2019). In recent literature, the performance of MM-PBSA and MM-GBSA rescoring (E. Wang et al., 2019; Z. Wang et al., 2019; Weng et al., 2019) have been evaluated also in protein–protein docking showing that MM-GBSA may be a good choice for predicting the binding affinities and identifying correct binding structures (Chen et al., 2016). Within this context, we’ve shown how the MM-GBSA rescoring is able to further improve the accuracy of the docking predictions. Recently, the same procedure of MM-GBSA rescoring has been applied to the three top scored docking poses showing that the results have improved compared the case in which only the best scored docking pose is considered (Rastelli et al., 2010; Sun, Li, Shen, et al., 2014). The method is a good compromise between efficiency and speed since it has been applied on minimized protein–ligand complexes. Future challenges could consist in finding more rigorous rescoring methods in order to improve the accuracy with a reasonable computational effort. Within this framework, the integration of binding pose prediction and affinity estimations based on artificial intelligence techniques could be implemented in order to improve the accuracy and sensitivity of the presented blind docking protocol.

5. Conclusions

The docking protocol developed in this study has demonstrated high efficiency in predicting protein–ligand docking when binding sites are a priori not known. It is worth to mention that the FRAD approach is a simple, fast, and easy to implement. FRAD is accurate in term of binding site identification and RMSD of the lowest energy docked pose. Nevertheless, the idea of performing a re-scoring on the FRAD docking results by methods estimating the free binding energy such as MM/GBSA proved to be a good solution for improvement. The adoption of the FRAD protocol has brought consistent increase of the ligand pose prediction accuracy with a reasonable computational effort, bringing an overall improvement in comparison to the BD procedure.

Authors’ contributions

GG and MAD conceived the research. GG, and ADG did the molecular dynamics simulations. GG, ADG, BM, DP, AD and MAD analyzed and rationalized the data. All authors wrote the paper and critically commented to the manuscript. All authors read and approved the final manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Funding

The present study has been developed as part of the VIRTUOUS project funded by the European Union’s Horizon 2020 research and innovation program under the Maria Sklodowska-Curie-RISE Grant Agreement No 872181 (https://www.virtuoush2020.com/). This work was also supported by a grant from the Swiss National Supercomputing Centre (CSCS).

Additional information

Funding

The present study has been developed as part of the VIRTUOUS project funded by the European Union’s Horizon 2020 research and innovation program under the Maria Sklodowska-Curie-RISE Grant Agreement No 872181 (https://www.virtuoush2020.com/). This work was also supported by a grant from the Swiss National Supercomputing Centre (CSCS).