Comparing the Performance of Molecular Docking Tools for HIV-1 Protease Inhibitors

. Human immunodeficiency virus (HIV) continues to pose a significant public health threat worldwide, disproportionately affecting marginalized and vulnerable populations despite advancements in prevention and treatment. One promising approach in the development of antiretroviral therapies is the inhibition of HIV-1 protease, a key enzyme in the virus life cycle. To this end, we evaluated 23 known inhibitors of HIV-1 protease using three docking methods - Autodock 4.2 (AD4), AutoDockVina (Vina), and a modified Vina (mVina)- and found that all three methods produced reasonable binding affinities for these ligands, but Vina performed best in terms of precision and docking success rates. The results provide important guidance for the selection of appropriate support tools for screening potential HIV-1 protease inhibitors in treating HIV/AIDS.


Introduction
Human Immunodeficiency Virus (HIV) is a viral infection that attacks the body's immune system, gradually weakening it over time. The virus is primarily transmitted through sexual contact, sharing of contaminated needles, and from mother to child during pregnancy, childbirth, or breastfeeding. If left untreated, HIV can progress to acquired immunodeficiency syndrome (AIDS), which is the most advanced stage of the infection. AIDS occurs when the immune system is severely damaged, making it vulnerable to opportunistic infections and cancers (Simon, Ho, & Karim, 2006). The virus has had a significant impact on global health, with more than 38 million people living with HIV worldwide (Montaner et al., 2006). Although there have been major strides in the development of antiretroviral therapy to manage the virus, HIV/AIDS remains a critical public health issue that requires continued attention and efforts towards prevention and treatment. HIV is a complex retrovirus that has a genome composed of two copies of a single-stranded RNA molecule (Deeks, Overbaugh, Phillips, & Buchbinder, 2015). The viral genome encodes for several key proteins, including reverse transcriptase, integrase, and protease, which are essential for the replication and survival of the virus. HIV is known for its high genetic variability, which results from a combination of high mutation rates during replication and recombination between different viral strains. This genetic diversity has led to the emergence of several distinct subtypes and recombinant forms of HIV, each with their own unique biological characteristics and geographic distribution (Santoro & Perno, 2013). Understanding the genomic structure and evolution of HIV is critical for the development of effective diagnostic tools, therapeutics, and vaccines to combat this global health threat. HIV-1 protease is a key enzyme in the life cycle of HIV, playing a crucial role in the maturation of new viral particles (Brik & Wong, 2003). The enzyme is responsible for cleaving the precursor proteins produced by the virus into their functional components, which are required for viral replication and infectivity. Inhibition of HIV-1 protease has been a successful strategy in the development of antiretroviral therapies for the treatment of HIV infection (Debouck, 1992). Several classes of HIV-1 protease inhibitors (PIs) have been developed, including peptidomimetic inhibitors and non-peptidic inhibitors, which target the active site of the enzyme and prevent its function (Deeks, Smith, Holodniy, & Kahn, 1997;Ghosh, Osswald, & Prato, 2016;Vacca & Condra, 1997;Wlodawer & Erickson, 1993). Although PIs have been shown to be highly effective in suppressing viral replication and reducing the risk of disease progression, the emergence of drug-resistant strains has highlighted the need for continued research into alternative therapeutic strategies (Ohtaka & Freire, 2005;Ridky & Leis, 1995). Nevertheless, HIV-1 protease remains an important and validated target for the development of new therapies for the treatment of HIV infection. In silico screening using computational approaches is a powerful tool for the discovery of new inhibitors. This approach enables the efficient screening of large compound libraries for potential inhibitors, reducing the time and cost associated with traditional experimental screening methods. Several computational methods have been developed for the identification of new HIV-1 protease inhibitors (PIs) (Fong, McNamara, Hillier, & Bryce, 2009;Han, Pelletier, & Hodge, 1998;Holloway et al., 1995;Ibrahim, A Saleh, M Elshemey, & A Elsayed, 2012;Judd et al., 2001;Lunney et al., 1994;Surleraux et al., 2005;S. Wang et al., 1996). These methods have been successfully used to identify new leads for the development of HIV-1 PIs with improved potency, selectivity, and resistance profiles. The integration of computational and experimental methods is an effective strategy for the development of new HIV-1 PIs for the treatment of HIV infection. Molecular docking is a computational technique used to evaluate the binding affinity of a large number of tested ligands to a target protein (Irwin & Shoichet, 2016;McDonnell, Reynolds, & Fogel, 1995;Sousa, Fernandes, & Ramos, 2006;Zhang, Li, Yu, & Jin, 2022). However, a critical aspect of any docking application is the scoring function, which determines the likelihood of a particular ligand-protein complex forming a stable structure. There are numerous servers, software packages, and applications available that use different force fields and algorithms to evaluate the scoring function.
Among the free open-source docking tools available, Autodock4 (AD4) (Morris et al., 2009), Autodock Vina (Vina) (Trott & Olson, 2010), and mVina (an experimental parameter correction version of Vina) (Pham et al., 2022) are notable for their ability to determine the ligand binding affinity quickly and accurately. These tools have been widely used in various drug discovery projects and have been shown to provide reliable results. In this study, we focused on comparing the performance of three docking software packages including AD4, Vina, and mVina. Specifically, we evaluated their efficacy in docking known HIV-1 PIs, with the aim of identifying an optimal approach for discovering new drugs to block HIV infection. Our findings are critical in selecting the most appropriate docking methodology for future drug discovery efforts in this area, potentially leading to the development of novel and effective AIDS treatments.

Structural Analysis
The GROMACS tools (Abraham et al., 2015) were used to calculate the root-mean-square deviation (RMSD) between two structures. The interactions of the docked poses of PI-HIV-1 protease complexes were visually represented using Pymol software.
Tab. 1: Free binding energy predicted by three docking software packages. The experimental affinity ∆Gexp was estimated via the equation ∆G = RT lnk i . where k i is the inhibition constant, R is the gas constant and T is the absolute temperature.

Estimated Ligand-Binding Free Energy
In Table 1, we present the free binding energy calculations for 23 ligands with HIV-1 PI using three software packages. We considered only the results with the lowest free binding energy or highest binding affinity due to their higher probability of actual occurrence. Overall, the three software packages produced reasonable docking scores within the range of -21.3 to -7.77 kcal.mol −1 , with most of the free binding energies being below -9.0 kcal.mol −1 . The results reveal that AD4 had the highest correlation coefficient (R =-0.47) between the experimental affinity ∆G exp and the predicted free binding energy, followed by mVina (-0.34) and Vina (-0.30), respectively. However, it's important to note that these three methods did not account for complex dynamics, explicit solvent effects, or the limited docking position of the ligand, leading to a small correlation in the prediction of ligand binding affinity.
When evaluating the accuracy of a method, the root-mean-square error (RMSE) is a key metric. A smaller RMSE suggests a higher level of accuracy. Table 2 presents the RMSE values for the three docking methods used when comparing the experimental and docking data. The results demonstrate that Vina outperformed the other two methods in terms of precision, while mVina exhibited the lowest accuracy. Furthermore, the standard error (SE) of the mean binding-free energy ∆G exp was calculated for each method, and the findings indicate that Vina produced the most reliable results, with the smallest SE, followed by mVina and AD4, respectively. Thus, the Vina docking method demonstrated faster convergence compared to the other two methods, further highlighting its efficacy and accuracy.

Successful docking rate
The success rate of docking was determined based on the similarity between the experimental ligand structure and the docking conformation, where a binding conformation was considered successful if the RMSD of atomic positions compared to the corresponding experimental structure was less than 0.2 nm (Bursulaya, Totrov, Abagyan, & Brooks, 2003). Table  3 shows the results of calculating the RMSD values, which indicate that Vina had the smallest mean RMSD of 0.13 ± 0.03 nm and that mVina had a similar average RMSD of 0.16 ± 0.03 nm. These results suggest that both Vina and mVina were able to accurately generate binding poses for HIV-1 PIs. In contrast, the mean RMSD value between the experimental structures and the docked poses using the AD4 protocol was 0.27 ± 0.03 nm, indicating that AD4 did not produce proper binding poses for the ligands. Note that the reported error is the standard error of the mean.
The suitability of a ligand-binding pose identification method is directly proportional to its docking success rate, as it represents the likelihood of finding the native binding pose. The success rate varies depending on the RMSD values, with smaller values indicating greater accuracy. Vina outperformed all other tested methods with a success rate of up to 83% and 59% at RMSD values below 0.15 and 0.1 nm, respectively ( Figure 1). mVina had success rates of 74% at RMSD values below 0.2 and 0.15 nm. However, the docking success rate of mVina fell to 39% when the RMSD was less than 0.1 nm. AD4 had poor success rates of only 30% and 9% at RMSDs below 0.2 and 0.1 nm, respectively, indicating that mostly it cannot accurately generate the binding pose of HIV-1 PIs. Overall, Vina consistently had a success rate above 50%, making it the superior method among the three approaches.

Redocked binding pose
Commonly used to validate docking results, molecular dynamics (MD) simulations are closely related to ligand-binding position accu-racy, a crucial factor in studying ligand-receptor interactions. The accuracy of MD simulations' free energy calculations depends on the input molecular structure. If the docking position differs significantly from the initial position, the MD simulation may need longer to reach equilibrium, and inadequate system control may lead to a false equilibrium. This greatly reduces binding affinity estimation accuracy. Therefore, evaluating the docking position, ranking the docking poses against experimental data, and determining which software generates docking poses closest to the original is necessary. This enables subsequent MD simulations to be more accurate and shorter. Vina and mVina were shown to produce suitable docking binding poses for HIV-1 PIs. Therefore, the binding structures of representative HIV-1 PIs gained by the two methods are subsequently analyzed. A high match be- tween the ligands from X-ray crystal structure and docking simulations expresses compatibility between the two positions ( Figure 2). The findings further strengthen the point that Vina is the most suitable docking tool for HIV-1 PIs among the three tested packages.
The experimental ligand pose (in green) was compared to the docking Vina (in wheat) or mVina (in light pink) poses.

Discussions
The question of which docking program is best is a common query in the field. In order to guide the selection of a suitable program for a specific application, several studies have been conducted to compare the benefits and drawbacks of various software, using different benchmarks as general reference points (Cross et al., 2009;Kellenberger, Rodrigo, Muller, & Rognan, 2004;Z. Wang et al., 2016). Both AD4 and Vina are frequently used docking programs that offer Windows compatibility. AD4 has been instrumental in discovering several potent inhibitors for peptides, proteins, and genes, while Vina has gained acceptance due to its user-friendliness and reliability in determining ligand-binding affinity (Gaillard, 2018). Vina's impressive computing capabilities have enabled it to predict the binding pose of large substrates to protein targets (Caffalette, Corey, Sansom, Stansfeld, & Zimmer, 2019;Vu et al., 2019) and estimate the binding affinities of small compounds to biomolecular targets (Grither & Longmore, 2018;Noike et al., 2015). However, Vina's  high docking success rate does not always correspond to a strong correlation between predicted and experimental binding affinity (Nguyen et al., 2020), making it challenging to rank ligandbinding affinity accurately. In response, mVina was developed with experimental parameter correction to enhance the program's ability to rank binding affinity (Pham et al., 2022).
Docking calculations have undoubtedly contributed to significant cost and time reductions in new drug research. However, like any other technique, molecular docking has its limitations. It is crucial to keep in mind that, although there are several reliable docking programs available, not all docking algorithms are suitable for every system. Proper selection of a suitable docking program based on the system's specific requirements is critical to obtaining accurate and reliable results (Cole, Davis, Jones, & Sage, 2017). Previously, docking experiments used AD4 and Pardock programs to study 25 HIV-1 protease-inhibitor complexes revealed that Pardock exhibited a higher correlation coefficient of 0.801 for predicted binding energy, whereas AD4 achieved a lower coefficient of 0.484 (Gujjula, 2008). Recently, the performance of Vina, 1-Click Docking, and UCSF DOCK was evaluated using 8 receptor-ligand complexes of HIV-1 protease inhibitors, indicating that Vina possessed the highest Pearson correlation coefficient of 0.81 with experimental results (Salih, 2022).
In our investigation, we re-docked 23 HIV-1 protease inhibitors using AD4, Vina, and mVina. Our analysis showed that for the particular HIV-1 protease inhibitors tested, all three methods produced reasonable ligand-binding affinities, but Vina demonstrated superior precision and docking success rates compared to the other two. Moreover, Vina was able to generate a ligand-binding pose that closely resembled the crystal structure obtained from practical experiments. It is recommended to conduct future research with a larger number of inhibitors, if they are available, to further support the results. The findings are then significant in guiding the selection of an appropriate support tool for screening potential HIV-1 protease inhibitors, which play a critical role in treating HIV/AIDS. Gujjula, K. R. (2008). Prediction and comparison of HIV-1 protease inhibitor binding energies by various molecular docking methods.