Protein–ligand affinity prediction via Jensen-Shannon divergence of molecular dynamics simulation trajectories

Abstract

Predicting the binding affinity between proteins and ligands is a critical task in drug discovery. Although various computational methods have been proposed to estimate ligand target affinity, the method of Yasuda et al. (2022) ranks affinities based on the dynamic behavior obtained from molecular dynamics (MD) simulations without requiring structural similarity among ligand substituents. Thus, its applicability is broader than that of relative binding free energy calculations. However, their approach suffers from high computational costs due to the extensive simulation time and the deep learning computations needed for each ligand pair. Moreover, in the absence of experimental ΔG values (oracle), the sign of the correlation can be misinterpreted. In this study, we present an alternative approach inspired by Yasuda et al.’s method, offering an alternative perspective by replacing the distance metric and reducing computational cost. Our contributions are threefold: (1) By introducing the Jensen–Shannon (JS) divergence, we eliminate the need for deep learning-based similarity estimation, thereby significantly reducing computation time; (2) We demonstrate that production run simulation times can be halved while maintaining comparable accuracy; and (3) We propose a method to predict the sign of the correlation between the first principal component (PC1) and ΔG by using coarse ΔG estimations obtained via AutoDock Vina.

This study proposes a method to estimate protein-ligand binding affinity using molecular dynamics (MD) simulation trajectories. Structural fluctuations of protein-ligand complexes are quantified by computing the Jensen-Shannon (JS) divergence between trajectory-derived distributions. The resulting similarity matrix is subjected to principal component analysis (PCA), where the primary components reflect the diversity in protein responses to different ligands. By projecting each ligand onto this reduced space, we investigate correlations between the PCA axes and experimental binding free energies. This framework enables rapid comparison of ligand-induced protein dynamics and offers a computationally efficient approach to prioritize compounds based on their dynamic impact.

Significance

This work proposes a complementary approach to existing MD simulation-based methods for binding affinity prediction, reducing computational costs and providing a new perspective on trajectory similarity through JS divergence. The proposed use of JS divergence and docking-based estimation of principal component polarity enhances both efficiency and robustness in scenarios where experimental binding data are unavailable.

 Introduction

Protein–ligand binding affinity prediction is a cornerstone technique in drug design and biological research. Binding affinity information is essential for screening candidate compounds that either bind to a specific protein or target a specific ligand. However, large-scale screening using biochemical experiments or structural analyses requires significant time and effort. Thus, computational approaches for rapid and efficient screening are highly desirable [1].

Existing computational methods for binding affinity prediction can be broadly classified into physics-based approaches and data-driven approaches. Physics-based methods, such as absolute free energy calculations [2], free energy perturbation (FEP) methods [3,4], and MM/PBSA [5], can yield highly accurate predictions but at a high computational cost. In contrast, data-driven approaches using machine learning or deep learning techniques [6–9] are typically faster during inference but require large, high-quality training datasets and are often subject to biases present in the training data.

An alternative strategy evaluates binding affinity by comparing the dynamic behavior of proteins upon ligand binding using molecular dynamics (MD) simulations [10]. In this approach, MD simulations are performed for both the apo form of the protein and its various ligand-bound complexes. The differences in the dynamics of binding site residues between ligand systems are quantified using Wasserstein distances computed via deep learning [11]. The resulting distance matrix is embedded into a three-dimensional space using nonlinear dimensionality reduction, followed by principal component analysis (PCA). Finally, the value along the first principal component (PC1) is used to predict lower ΔG values. Yasuda et al. demonstrated good correlations between PC1 and experimental ΔG for targets such as Bromodomain 4 (BRD4) (R=0.88) and protein phosphatase 1B (PTP1B) (R=0.70). Nonetheless, the method is hampered by high computational costs from both the MD simulations and the deep learning computations required for each ligand pair. Moreover, even if PC1 values can be obtained, the sign of the correlation with ΔG remains ambiguous without oracle data.

In this study, we address these issues by reducing computational cost and expanding the method’s applicability. Specifically, we avoid the deep learning-based Wasserstein distance computation by employing the Jensen–Shannon (JS) divergence for estimating similarity between trajectories. In addition, we propose a method to determine the sign of the correlation between PC1 and ΔG using coarse ΔG estimations from docking simulations.

While the approach of Yasuda et al. demonstrated promising correlations between PC1 and experimental ΔG, it requires computationally intensive procedures, such as long MD simulations and deep learning-based trajectory comparisons. Moreover, since the polarity of the PC1–ΔG correlation depends on prior knowledge of binding affinities, the interpretation may be ambiguous in practice. Our method addresses these points by simplifying the similarity metric and introducing a practical means to estimate the correlation direction without relying on exhaustive experimental data.

 Materials and methods

Figure 1 presents an overview of the proposed method. While the overall framework is similar to the method of Yasuda et al. [10], the key differences are in the computation of the distance matrix (step d) and the estimation of the sign of the PC1 correlation (step g). In our method, deep learning is replaced by the JS divergence for rapid evaluation of similarity between trajectories.

Figure 1  Overview of the proposed method. The main differences compared to the method of Yasuda et al. [10] are the computation of the distance matrix W (step d) and the estimation of the polarity of PC1 (step g).

 Dataset

We evaluated the proposed method using three datasets: BRD4 and PTP1B as used in previous studies, and c-Jun N-terminal kinase (JNK1), which has been employed as a benchmark in FEP calculations [12]. These datasets are widely used in free energy calculations [2,4] and are appropriate for evaluating the proposed method. The crystal structures of the complexes and the chemical structures of the compounds are provided in Supplementary Figures S1–S4.

In this study, we refer to experimentally determined binding affinities as the “oracle ΔG,” following a convention used in computational benchmarking studies to denote ground-truth or reference values. These values serve as the standard for evaluating the predictive performance of our method. Experimental binding free energy (ΔG) values used for evaluation in this study were obtained from previously published literature. Specifically, BRD4 affinities were taken from Aldeghi et al. [2], PTP1B affinities from Wilson et al. [13], and JNK1 affinities from Szczepankiewicz et al. [12]. In contrast, ΔG_dock refers to the docking score obtained from AutoDock Vina [14] for the top-ranked pose of each ligand. This score is not an exact free energy value but a heuristic estimate reflecting the predicted binding strength based on a scoring function that includes empirical and knowledge-based energy terms. All docking simulations were performed using default AutoDock Vina parameters.

 Preparation of initial protein–ligand complex structures

We used a total of 11 ligands for BRD4, 8 for PTP1B, and 10 for JNK1 in this study (see Supplementary Figures S1–S4 for chemical structures).

For BRD4, the initial complex structures were prepared according to the procedure described by Aldeghi et al. [2]. Specifically, the initial conformations were taken from holo crystal structures (PDB IDs: 3U5J, 3U5L, 4OGI, 4OGJ, 3MXF, 4MR3, 4MR4, 3SVG, 4J0R, 4HBV) and from docking results of the ligand into the apo protein (PDB ID: 2OSS). For the docked complexes, the docking was performed so that the ligand-binding moieties overlapped with those observed in the reference holo structures. All complex structures were aligned to ensure consistency prior to simulation.

For PTP1B and JNK1, the initial ligand-bound conformations were prepared based on the protocol described in Wang et al. [4]. PTP1B used the apo crystal structure with PDB ID 2QBS, and JNK1 used the apo structure with PDB ID 2GMX. The ligand positions for both systems were determined by superimposing the ligand structures onto a reference complex, ensuring consistent orientation and placement within the binding pocket.

Hydrogen atoms were added to each PDB file using H++ [15] at pH 7.0 to complete missing atoms.

 Molecular dynamics simulation

All-atom MD simulations were conducted for both the protein and the protein–ligand complexes. Each hydrated system was subjected to energy minimization, followed by restrained MD simulations in the NVT and NPT ensembles, and finally a production run.

A cubic box with a 10.0 Å buffer region was created using the TIP3PBOX model [16] for solvation. The ff14SB and GAFF force fields were applied for the protein and ligand, respectively. The system charge was neutralized by adding sodium or chloride ions, and topology and coordinate files were generated.

Using Amber22 [17], each system was first minimized using 5,000 steps of steepest descent. Subsequently, a 100 ps NVT simulation at 300 K was performed with a restraint of 10 kcal/mol·Å² on heavy atoms. This was followed by a 100 ps NPT simulation under the same restraint at 300 K and 1 bar. Finally, a 400 ns production run was performed under the NPT ensemble with random initial velocities, and trajectories were saved every 2 ps. Each system was simulated in three independent trials.

 Identification of binding site residues

Binding site residues were identified based on the activity ratio, defined as n/N, where n is the number of frames in which the minimum distance between any heavy atom of an amino acid residue and any heavy atom of the ligand is less than 5 Å, and N is the total number of frames. Residues with an activity ratio greater than 0.5 during the first 100 ns of simulation were considered part of the binding site, and the union of such residues across all ligand complexes was extracted (see Supplementary Figure S5). The number of identified binding site residues was 13 for BRD4, 17 for PTP1B, and 18 for JNK1. These residues were used for all subsequent trajectory and similarity analyses. The RMSD trajectories for binding site residues in each production run are provided in Supplementary Figures S6, S7, and S8. From the trajectories of binding site residues, rotation and translation were removed by fitting the trajectories to an identical structure in the backbone atoms of the binding site residues.

 Trajectory distance calculation

For each protein or protein–ligand complex, the trajectories of the binding site residues were used to estimate probability density functions via kernel density estimation using a Gaussian kernel. The bandwidth was determined using Silverman’s rule-of-thumb. The similarity between probability density functions from different ligand systems was then evaluated using the Jensen–Shannon (JS) divergence.

For systems i and j, the distance W_ij is defined as

  

$$\begin{array}{cc}\hfill W_{ij}\hfill & \hfill =\frac{1}{N_a}\sum _{a=1}^{N_a}J{S}_a(i\parallel j),\hfill \\ \hfill JS(i\parallel j)\hfill & \hfill =\frac{1}{2}(KL(i\parallel M)+KL(j\parallel M)),\hfill \\ \hfill M(𝐱)\hfill & \hfill =\frac{1}{2}(i(𝐱)+j(𝐱)),\hfill \\ \hfill KL(i\parallel j)\hfill & \hfill =\sum _{𝐱\in 𝒢}^{}i(𝐱)log\frac{i(𝐱)}{j(𝐱)},\hfill \end{array}$$

where KL(· || ·) denotes the Kullback–Leibler divergence. Here, i(x) and j(x) represent the probability density functions estimated via kernel density estimation from the trajectories of ligand systems i and j, respectively. The variable x denotes the spatial coordinate (x, y, z) corresponding to the centroid of heavy atoms of each binding site residue.

The summation is taken over a discrete spatial grid $𝒢$, which is defined as a regular 3D mesh covering the union of binding site coordinates across all systems. The grid spacing along each axis (x, y, and z) was set to 1/100 of the range of aligned centroid coordinates along that axis, computed over all frames and systems. The number of binding site residues N_a depends on the protein system; for example, in BRD4, N_a=13.

 Binding affinity evaluation based on dynamic similarity

The distance matrix W={W_ij} was embedded into a low-dimensional vector space. A set of 5-dimensional vectors was determined by minimizing the error function

  

$${𝐩}_1,{𝐩}_2,\dots ,{𝐩}_L=argmin\sum _{i<j}^{}(W_{ij}-\parallel {𝐩}_i-{𝐩}_j\parallel ),$$

using simulated annealing and quasi-Newton methods. Here, L denotes the dimensionality of the embedded vectors, which was set to 5 in this study. PCA was then applied to the embedded vectors, and the correlation between the first principal component (PC1) and the binding free energy ΔG was evaluated.

 Estimation of the sign of the PC1–ΔG correlation

In the previous method [10], the correct orientation of PC1 with respect to ΔG was ambiguous without knowledge of the oracle ΔG. In our approach, a pseudo-oracle is obtained via docking simulations using AutoDock Vina [14]. For each protein–ligand pair, the docking score is calculated and its correlation with PC1 is determined. If PC1 is positively correlated with the docking score, we assume a positive correlation with ΔG, and vice versa. Although the docking score does not precisely match the experimental ΔG, it is sufficiently reliable to determine the polarity of the correlation.

 Computational environment

MD simulations were performed on the TSUBAME 4.0 GPU node (gpu 1 queue: AMD EPYC 9654 2.4 GHz [8 cores], NVIDIA H100 SXM5 94GB HBM2e [1 card], 96 GiB DDR5-4800 RAM) at the Center for Information Infrastructure, Institute of Science Tokyo. Amber22 [17] was used for MD simulations, and subsequent analyses were conducted in Python. The software libraries and versions used are listed below:

Amber	22 Update5
AmberTools	23 Update6
AutoDock Vina	1.2.5
Python	3.9.19
NumPy	1.26.4
MDAnalysis	2.7.0
scikit-learn	1.3.0
SciPy	1.12.0

 Results and discussion

 Comparison of complex dynamics using JS divergence

Figure 2 shows the distance matrix computed from the JS divergence between the probability distributions derived from the binding site trajectories. It is evident that the distances between the apo system (system 0) and the ligand-bound systems are larger compared to those among the complexes, indicating significant dynamic differences induced by ligand binding.

Figure 2  Distance matrix based on the JS divergence between probability distributions from the binding site trajectories. Higher (yellow) values indicate greater dynamic differences. System number 0 represents the apo protein trajectory.

 Principal component analysis of embedded vectors

The distance matrix W was embedded into a 5-dimensional vector space and subjected to PCA. As shown in Figure 3, for all datasets the ligand systems with high and low ΔG values are clearly separated along PC1.

Figure 3  PCA results after embedding the distance matrix W into a 5-dimensional space. Black dots indicate the apo (ligand-free) system.

 Binding affinity evaluation based on dynamic similarity

The correlation between PC1 and the oracle ΔG is presented in the upper panel of Figure 4. For the BRD4 and PTP1B datasets, PC1 exhibits a positive correlation with ΔG, whereas for the JNK1 dataset, a negative correlation is observed.

Figure 4  Correlation plots of PC1 with the oracle ΔG. The upper panel shows results using trajectories up to 400 ns, while the middle and lower panels correspond to shorter simulation durations.

Although we focused on the correlation between PC1 and ΔG following the approach of Yasuda et al. [10], it is possible that other principal components, such as PC2, also carry predictive signals. However, since the PCA space depends on the specific dataset and ligand distribution, incorporating multiple components into a unified regression model is not straightforward and may reduce generalizability. For this reason, we limited our analysis to PC1 in this study to maintain consistency and interpretability. Investigating multivariate models involving PC2 remains a subject for future work.

In the original study by Yasuda et al. [10], both examined systems exhibited the same direction of correlation between PC1 and ΔG, which supported a consistent interpretation of PC1 as an affinity-associated axis. However, in our results, the JNK1 dataset showed an inverse correlation, indicating that the sign of the correlation is not necessarily preserved across different systems. This suggests that PC1 may not universally represent the same physical aspect of protein–ligand dynamics. At present, we do not have sufficient evidence to interpret the mechanistic cause of this inversion. Further stud-ies across a broader range of targets will be necessary to clarify whether PC1 consistently captures physically meaningful dynamic features linked to binding affinity.

 Effect of MD simulation length on prediction accuracy

We investigated the impact of production run simulation time by comparing trajectories of 0–400 ns, 0–200 ns, and 0–100 ns (Figure 4). For the BRD4 dataset, the separation of ligand systems by ΔG along PC1 is evident in all cases. However, for PTP1B and JNK1, while 200 ns and 400 ns trajectories yield clear separation between high and low ΔG systems, the 100 ns trajectories do not. These results suggest that a simulation time of around 200 ns or more is necessary to capture dynamic behavior reflecting structural stabilization and binding affinity differences. However, because the deep learning-based Wasserstein distance method used by Yasuda et al. is not publicly available, we were unable to reproduce their approach or examine whether similar reductions in simulation time would retain predictive accuracy. Nevertheless, the consistent separation observed in both studies indicates that trajectory-based methods may generally permit shortened simulations without substantial loss of information.

 Estimation of the PC1–ΔG correlation sign

Figure 5 (upper panel) compares PC1 with the docking scores ΔG_dock obtained via AutoDock Vina. For BRD4 and PTP1B, the signs of the correlation coefficients between PC1 and ΔG_dock match those between PC1 and the oracle ΔG, validating the use of docking scores for sign determination. For JNK1, however, an inverse sign was observed. The lower panel of Figure 5 shows that although the docking scores do not perfectly agree with the oracle ΔG, they are sufficiently reliable for determining the polarity of the PC1–ΔG correlation. In cases where a few compounds have been evaluated using absolute binding free energy calculations such as FEP or MP-CAFEE [18], those values could also be used to determine the sign.

Figure 5  (Upper panel) Relationship between PC1 and docking scores ΔG_dock. (Lower panel) Comparison between docking scores ΔG_dock and oracle ΔG.

 Comparison of prediction accuracy and computational cost

Table 1 compares the computational cost and performance between the method of Yasuda et al. [10] and the present work. In the previous approach, deep neural network-based Wasserstein distance calculations required several hours per ligand pair using an NVIDIA GeForce RTX 3090 24GB GPU, whereas our method employing JS divergence computes the distance in a matter of seconds. Although the correlation between PC1 and ΔG is comparable, our method reduces both the required simulation time (200 ns instead of 400 ns) and the per-comparison computational cost.

Table 1 Comparison with the previous study [10]

	BRD4		PTP1B
	previous method	this study	previous method	this study
correlation (PC1 vs. ΔG)	R=0.88	R=0.66	R=0.70	R=0.78
required simulation length	400 ns	200 ns	400 ns	200 ns
similarity computation	DNN-based Wasserstein distance	JS divergence	DNN-based Wasserstein distance	JS divergence
time per pair comparison	several hours	several seconds	several hours	several seconds

Additionally, we compared our method with relative FEP calculations for the PTP1B and JNK1 datasets [19]. Although direct comparison is challenging due to differences in simulation protocols, we estimated the total simulation times for the FEP calculations based on the number of λ windows and sampling lengths reported in a previous study [19]. Specifically, the number of λ windows was 656 for PTP1B and 504 for JNK1, with each window involving approximately 500 ps of equilibration and 4 ns of production run using a 4.0 fs time step. This results in estimated total simulation times of approximately 3,000 ns for PTP1B and 2,300 ns for JNK1, respectively. In comparison, the total simulation time in our study was approximately 6,800 ns for PTP1B and 7,600 ns for JNK1. Despite the protocol differences, the correlation (R²) values were comparable (FEP: R²=0.47 for PTP1B and R²=0.74 for JNK1; this study: R²=0.61 for PTP1B and R²=0.48 for JNK1). It should be noted that while relative FEP directly computes free energy differences (ΔΔG), our method provides ranking based on overall dynamic similarity without requiring structural overlap or perturbation pathways.

 Conclusion

In this study, we have presented a modified framework for protein–ligand binding affinity prediction, inspired by Yasuda et al.’s method but based on a distinct similarity metric and analysis approach. By introducing the Jensen–Shannon divergence, we can rapidly compute the distance matrix between system trajectories without the need for deep learning, thereby significantly reducing computational cost. We also demonstrated that the production run simulation time can be reduced from 400 ns to 200 ns while maintaining comparable accuracy. Finally, by incorporating coarse ΔG estimations from AutoDock Vina, we proposed a strategy to determine the sign of the PC1–ΔG correlation.

Although our method requires three-dimensional structural information of the target, the increasing accuracy of structure prediction tools such as AlphaFold3 [20] suggests that our approach may soon be applicable even to targets with unresolved experimental structures. Future work will explore the application of our method using predicted structures.

 Conflict of interest

The authors declare that they have no conflict of interest.

 Author contributions

K.I. was responsible for executing the research, implementing and analyzing the work, visualizing the results, drafting the article, and revising the manuscript. M.O. provided overall research supervision, conceived and designed the work, led the discussion of the results, contributed to drafting the article, and revised the manuscript. All authors approved the final version to be published.

 Data availability

The datasets, trajectories, and estimation code (python) can be obtained from the following link; https://drive.google.com/drive/folders/1mxiErt0jSKVACIVkahntLoS82tZFDow7/. A preliminary version of this work was deposited in bioRxiv (https://doi.org/10.1101/2025.02.17.638772) on February 23, 2025.

 Acknowledgements

This work was financially supported by the Japan Science and Technology Agency FOREST (Grant No. JPMJFR216J), the Japan Society for the Promotion of Science KAKENHI (Grant Nos. JP23H04880, JP23H04887), and the Japan Agency for Medical Research and Development Basis for Supporting Innovative Drug Discovery and Life Science Research (Grant No. JP25ama121026). The authors thank Ikki Yasuda for valuable discussions regarding the computational cost of the deep neural network-based Wasserstein distance calculation, and Kairi Furui for contributions to the discussion of computational time in FEP calculations.

References

• [1]   Wang,  Y.,  Jiao,  Q.,  Wang,  J.,  Cai,  X.,  Zhao,  W.,  Cui,  X. Prediction of protein-ligand binding affinity with deep learning. Comput. Struct. Biotechnol. J. 21, 5796–5806 (2023). https://doi.org/10.1016/j.csbj.2023.11.009
• [2]   Aldeghi,  M.,  Heifetz,  A.,  Bodkin,  M. J.,  Knapp,  S.,  Biggin,  P. C. Accurate calculation of the absolute free energy of binding for drug molecules. Chem. Sci. 7, 207–218 (2016). https://doi.org/10.1039/C5SC02678D
• [3]   Cournia,  Z.,  Allen,  B.,  Sherman,  W. Relative binding free energy calculations in drug discovery: Recent advances and practical considerations. J. Chem. Inf. Model. 57, 2911–2937 (2017). https://doi.org/10.1021/acs.jcim.7b00564
• [4]   Wang,  L.,  Wu,  Y.,  Deng,  Y.,  Kim,  B.,  Pierce,  L.,  Krilov,  G., et al. Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field. J. Am. Chem. Soc. 137, 2695–2703 (2015). https://doi.org/10.1021/ja512751q
• [5]   Genheden,  S.,  Ryde,  U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin. Drug Discov. 10, 449–461 (2015). https://doi.org/10.1517/17460441.2015.1032936
• [6]   Gomes,  J.,  Ramsundar,  B.,  Feinberg,  E. N.,  Pande,  V. S. Atomic convolutional networks for predicting protein-ligand binding affinity. arXiv (2017). https://doi.org/10.48550/arXiv.1703.10603
• [7]   Jiménez,  J.,  Skalic,  M.,  Martinez-Rosell,  G.,  De Fabritiis,  G. K_DEEP: Protein–ligand absolute binding affinity prediction via 3D-convolutional neural networks. J. Chem. Inf. Model. 58, 287–296 (2018). https://doi.org/10.1021/acs.jcim.7b00650
• [8]   Öztürk,  H.,  Özgür,  A.,  Ozkirimli,  E. DeepDTA: Deep drug–target binding affinity prediction. Bioinformatics 34, i821–i829 (2018). https://doi.org/10.1093/bioinformatics/bty593
• [9]   Abbasi,  K.,  Razzaghi,  P.,  Poso,  A.,  Amanlou,  M.,  Ghasemi,  J. B.,  Masoudi-Nejad,  A. DeepCDA: Deep cross-domain compound–protein affinity prediction through LSTM and convolutional neural networks. Bioinformatics 36, 4633–4642 (2020). https://doi.org/10.1093/bioinformatics/btaa544
• [10]   Yasuda,  I.,  Endo,  K.,  Yamamoto,  E.,  Hirano,  Y.,  Yasuoka,  K. Differences in ligand-induced protein dynamics extracted from an unsupervised deep learning approach correlate with protein–ligand binding affinities. Commun. Biol. 5, 481 (2022). https://doi.org/10.1038/s42003-022-03416-7
• [11]   Endo,  K.,  Yuhara,  D.,  Tomobe,  K.,  Yasuoka,  K. Detection of molecular behavior that characterizes systems using a deep learning approach. Nanoscale 11, 10064–10071 (2019). https://doi.org/10.1039/C9NR00219G
• [12]   Szczepankiewicz,  B. G.,  Kosogof,  C.,  Nelson,  L. T.,  Liu,  G.,  Liu,  B.,  Zhao,  H., et al. Aminopyridine-based c-Jun N-terminal kinase inhibitors with cellular activity and minimal cross-kinase activity. J. Med. Chem. 49, 3563–3580 (2006). https://doi.org/10.1021/jm060199b
• [13]   Wilson,  D.,  Wan,  Z.,  Xu,  W.,  Kirincich,  S.,  Follows,  B.,  Joseph-McCarthy,  D., et al. Structure-based optimization of protein tyrosine phosphatase 1B inhibitors: From the active site to the second phosphotyrosine binding site. J. Med. Chem. 50, 4681–4698 (2007). https://doi.org/10.1021/jm0702478
• [14]   Eberhardt,  J.,  Santos-Martins,  D.,  Tillack,  A. F.,  Forli,  S. AutoDock Vina 1.2.0: New docking methods, expanded force field, and Python bindings. J. Chem. Inf. Model. 61, 3891–3898 (2021). https://doi.org/10.1021/acs.jcim.1c00203
• [15]   Gordon,  J. C.,  Myers,  J. B.,  Folta,  T.,  Shoja,  V.,  Heath,  L. S.,  Onufriev,  A. H++: A server for estimating pK_as and adding missing hydrogens to macromolecules. Nucleic Acids Res. 33, W368–W371 (2005). https://doi.org/10.1093/nar/gki464
• [16]   Jorgensen,  W. L.,  Chandrasekhar,  J.,  Madura,  J. D.,  Impey,  R. W.,  Klein,  M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935 (1983). https://doi.org/10.1063/1.445869
• [17]   Case,  D. A.,  Aktulga,  H. M.,  Belfon,  K.,  Ben-Shalom,  I. Y.,  Berryman,  J. T.,  Brozell,  S. R., et al. Amber 2022. University of California, San Francisco (2022).
• [18]   Okada,  O.,  Yamashita,  H.,  Takedomi,  K.,  Ono,  S.,  Sunada,  S.,  Kubodera,  H. Prediction of the binding affinity of compounds with diverse scaffolds by MP-CAFEE. Biophys. Chem. 180, 119–126 (2013). https://doi.org/10.1016/j.bpc.2013.07.005
• [19]   Furui,  K.,  Ohue,  M. Benchmarking HelixFold3-predicted holo structures for relative free energy perturbation calculations. ACS Omega 10, 11411–11420 (2025). https://doi.org/10.1021/acsomega.4c11413
• [20]   Abramson,  J.,  Adler,  J.,  Dunger,  J.,  Evans,  R.,  Green,  T.,  Pritzel,  A., et al. Accurate structure prediction of biomolecular interactions with AlphaFold 3. Nature 630, 493–500 (2024). https://doi.org/10.1038/s41586-024-07487-w