Precise binding free energy calculations for multiple molecules using an optimal measurement network of pairwise differences

In drug discovery research, candidate drug molecules are increasingly ranked and selected by their computationally predicted affinities for the biomolecular target of interest before experimental synthesis and testing ^{1, 2, 3, 4, 5, 6}. This allows for the exploration of many more molecular design ideas than is practical with experimental approaches alone. Alchemical binding free energy (BFE) calculations have become a valuable tool in predicting the protein-ligand binding affinities and their utility continues to increase thanks to advances in the methodologies ^{7, 8, 9, 10, 11, 12, 4, 13} force field ^{14, 15, 16, 17, 18, 19}, computational power afforded by the graphical prcessing units (GPUs) ^{11, 9, 8}, and the availability of simulation software and simplified workflows ^{11, 9, 4, 20, 21, 22, 23, 24, 25, 26}.

BFE calculations can either compute the binding free energy of an individual molecule, by a technique commonly referred to as "absolute" binding free energy (ABFE) calculations ^{27, 28, 29, 30}, or compute the difference between the binding free energies of two molecules, by "relative" binding free energy (RBFE) calculations ^{31, 32, 11, 33}. The binding free energies of a set of molecules can be determined by a combination of ABFE calculations of select molecules and RBFE calculations of select pairs of molecules ³⁴.

Statistical errors in BFE calculations confound decision-making in candidate molecule selection and the assessment of their predictive accuracy. There have been significant methodological progress in minimizing statistical errors in single BFE calculations ^{35, 36}, but works have only recently been published addressing the question of how to reduce the collective statistical errors for BFE calculations of a set of molecules ^{34, 37}, which represents the majority of applications of BFE calculations in drug discovery.

As suggested previously, at a fixed total computational cost, the overall statistical errors in BFE calculations for multiple molecules can be minimized by optimally allocating the computational resources to different RBFE and ABFE calculations. In particular, minimizing the total variance in the estimated binding free energies with respect to the allocation corresponds to the A-optimal in experimental design ³⁴. Such optimal allocations, however, require knowledge of the statistical fluctuations in individual BFE calculations, which are unknown a priori. Because of this reason, A-optimal BFE calculations have not yet been reported for realistic examples.

In this work, we develop a method based on the A-optimal design to improve the statistical precision of BFE calculations. Our method, Network Binding Free Energy (NetBFE), performs adaptive BFE calculations in multiple iterations. In each iteration, NetBFE estimates the statistical fluctuations in different BFE calculations based on the results from previous iterations, and it reallocates the computational resources accordingly. This allows the A-optimal allocation to be used in realistic BFE calculations. We also extend the previous theoretical work by deriving new equations for the A-optimal design and the covariance matrix when only RBFE calculations are included (the previous theory required the inclusion of at least one ABFE calculation), a common scenario in the lead optimization stage of drug discovery.

Here, we apply NetBFE to calculating the binding free energies of 16 ligands against the TYK2 protein ^{11, 38} and 8 ligands against the HIF-2$\alpha$ protein ³⁹. We compare the statistical errors in the estimated binding free energies from NetBFE calculations and those in the BFE calculations using RBFE pairs selected by the LOMAP algorithm developed by the Mobley lab ⁴⁰, which is widely used in such calculations. At the same total computational cost, NetBFE reduces the statistical variance by about half compared to LOMAP. Our results suggest that NetBFE can substantially improve the statistical precision of BFE calculations in drug discovery, which should accelerate the computational selection of candidate drug molecules and thereby increase the efficiency of the design-make-test cycle.

We first illustrate the iterative allocation in NetBFE consisting of both ABFE and RBFE calculations, using the small example of eight inhibitors against the protein kinase TYK2. No experimental binding free energies were included in the allocation and analysis, i.e. $\delta_{i}^{-1}=0$ in Eq. 4. The estimated fluctuations $\{s_{e}\}$ and the corresponding A-optimal allocation of $\{n_{e}\}$ in each iteration are shown in Fig. 1A and Fig. 1B, respectively. The total variance $\mathrm{tr}(\mathbf{C}(n^{(i)},s^{(i)}))$ decreases with increasing iterations.

How well does this iterative procedure approach the true A-optimal allocation? If we assume that the estimate of $\{s_{e}^{(i=5)}\}$ after the last iteration represents the true statistical fluctuations, $\mathrm{tr}(\mathbf{C}(n^{(5)},s^{(5)}))$ approximates the variance of the A-optimal allocation. We measure how far the allocation $n^{(i)}$ in the $i$’th iteration—normalized to $n^{\prime(i)}\equiv\sum_{e}n_{e}^{(5)}(\sum_{e}n_{e}^{(i)})^{-1}n^{(i)}$ so that its total number of samples is the same as that in the last iteration—is from the this approximate A-optimal by comparing $\mathrm{tr}(\mathbf{C}(n^{\prime(i)},s^{(5)}))$ to $\mathrm{tr}(\mathbf{C}((n^{(5)},s^{(5)}))$. The former indeed approaches the latter (Fig. 1C) after only three iterations, suggesting that the iterative allocation is converging quickly to the true A-optimal.

Because ABFE converges much slower than RBFE (on average ABFE has a much higher statistical fluctuation than RBFE: $\langle\{s_{i}\}\rangle_{\mathrm{ABFE}}/\langle\{s_{ij}\}\rangle_{\mathrm{RBFE}}$ $\approx$ 6.2), most of the computations (approximately 81$\%$ of total resources) are allocated to ABFE in order to improve the precision of the absolute values of the computed binding free energies.

Often it is unnecessary to predict the absolute values of the binding free energies, such as in the case when the values for some similar molecules have been experimentally determined, thus the values for the molecules of interest can be determined by predicting the differences between the molecules, or in the case of ranking the molecules by their binding affinities, such that only their relative differences will suffice. In such cases ABFE calculations may be omitted to avoid the most expensive part of the computations, as noted above. In the following, we consider NetBFE with only RBFE calculations.

We used NetBFE to compute the binding free energies of 16 inhibitors against TYK2 ^{38, 11} and 8 inhibitors against HIF-2$\alpha$ ³⁹. These molecules are depicted in Fig. 2 and listed in the Supplementary Materials. We first characterized the progression of statistical errors after each iteration of extending the total simulation time (Fig. 3). The standard deviations and the total variance by and large decreases with increasing total simulation time. Moreover, the distribution of the standard errors becomes narrower, indicating that most of the free energy results have converged to similar precision.

In the case of HIF-2$\alpha$ inhibitors, the agreement between the experimental binding free energies and their predicted values improves at each iteration (Fig. 3), indicating that the improvement in precision is accompanied by the improvement in accuracy; for TYK2 inhibitors, the accuracy is quite good (RMSE < 1 kcal/mol) in all 5 iterations. Of course, accuracy in the binding free energy calculations depends on other factors such as the force field parametrization and the consistency between the computational setup and the experimental conditions. It remains unknown whether in other protein-ligand combinations or in other force fields the accuracy of the predicted binding free energies improves with the statistical precision.

In each iteration of the NetBFE calculation, the allocations $\{n_{e}\}$ are optimized according to the estimate of statistical fluctuations $\{s_{e}\}$ based on the results of the previous iteration, thus both $\{s_{e}\}$ and $\{n_{e}\}$ change from one iteration to the next. We measure the difference between two allocations of $\{n_{e}\}$ and $\{n^{\ast}_{e}\}$ by the Kullback-Leibler (KL) divergence:

where $f_{e}\equiv n_{e}/\sum_{e^{\prime}}n_{e^{\prime}}$ and similarly $f^{\ast}_{e}$ are the fractional allocations. In the two cases above, $\{n_{e}\}$ approach their limiting values with the progression of iterations (Fig. 4): the KL-divergence between $\{n_{e}^{(i)}\}$ and $\{n_{e}^{(\mathrm{opt})}\}$ decreases toward 0.69 for TYK2 and 0.28 for HIF-2$\alpha$ with each iteration, where $\{n_{e}^{(\mathrm{opt})}\}$ is the the A-optimal allocation based on $s^{(i=5)}$ from the last iteration, assuming, as in the above, that the estimate of $s^{(i=5)}$ from the last iteration is a good approximation to the true statistical fluctuations. Taking the A-optimal variance to be $\mathbf{C}_{\mathrm{opt}}\approx\mathbf{C}(n^{(\mathrm{opt})},s^{(5)})$, $\mathrm{tr}(\mathbf{C}(n^{\prime(i)},s^{(5)})/\mathrm{tr}(\mathbf{C}_{\mathrm{opt}})$ rapidly approaches one with each iteration for both TYK2 and HIF-2$\alpha$ examples (as above, $n^{\prime(i)}=\sum_{e}n_{e}^{(5)}(\sum_{e}n_{e}^{(i)})^{-1}n^{(i)}$ is the normalized allocation in iteration $i$). Note that $n^{(i)}$ is determined on a 2-connected sub-network (see Section 2.4) while $n^{(\mathrm{opt})}$ is determined on the full network. These observations suggest that the iterations in NetBFE using a two-connected sub-network result in near-optimal allocations.

A simple consistency check can help identify problematic RBFE calculations by comparing the $\Delta\Delta G_{ij,\mathrm{RBFE}}$ in the RBFE calculation against the difference $\Delta G_{i,\mathrm{NetBFE}}-\Delta G_{j,\mathrm{NetBFE}}$, where $\Delta G_{i,\mathrm{NetBFE}}$ and $\Delta G_{j,\mathrm{NetBFE}}$ are derived by the maximum likelihood estimator of Eq. 3. Fig. 5 shows such comparison for the NetBFE calculations of TYK2 and HIF-2$\alpha$ inhibitors; points significantly deviating from the diagonal (e.g. points for which $|\Delta\Delta G_{ij,\mathrm{RBFE}}-(\Delta G_{i,\mathrm{NetBFE}}-\Delta G_{j,\mathrm{NetBFE}})|\geq 3\sqrt{\sigma_{ij}^{2}+\sigma_{i}^{2}+\sigma_{j}^{2}}$) indicate potentially problematic RBFE calculations.

Having demonstrated above that the iterative procedure in NetBFE allocates the computational resources approximately according to the A-optimal, we next investigate whether such allocations yield improved precision in the computed binding free energies. We performed $g=3$ independent replicates of NetBFE calculations of the TYK2 and HIF-2$\alpha$ examples, initializing each replicate with the same binding poses but different random number seeds in equilibration and alchemical simulations, and computed the statistical variances in the predicted binding free energies across the replicates. In comparison, we computed the same sets of binding free energies—and their statistical variances from $g=3$ independent replicates—by allocating the same total computational cost evenly to pairs of RBFE calculations selected by a previously published method, LOMAP ⁴⁰. Fig. 6 shows the distributions of variances for the two NetBFE allocation methods (NetBFE($N$) and NetBFE($\Delta N$), explained in Section 2.4) and the LOMAP allocation method. The mean and median values of the variance distributions are listed in Table 1. The NetBFE allocations indeed significantly outperform the LOMAP allocations in terms of the statistical precision in the predicted binding free energies. The reduction in the statistical variance by approximately a factor of two is consistent with the previous results from simulated data ³⁴.

Reference binding free energy values from experimental measurements for a subset of the molecules, when available, may be incorporated in the NetBFE calculations by supplying $\tilde{x}_{i}$ and $\delta_{i}^{-1}>0$ in Eq. 4, 5. Here we tested the allocation and the corresponding precision and accuracy in NetBFE calculations that include different numbers of reference values. We performed four independent NetBFE calculations for the 16 inhibitors against TYK2, supplying reference experimental values for $k=0,1,2,4$ molecules, respectively.

The optimal allocations at the end of five iterations are shown in Fig. 7. When there are a small number of reference molecules ($k=1,2$), they serve as network hubs, in that there is a RBFE edge connecting them to every other molecule. As we increase the number of reference molecules ($k=4$), however, only a partial subset of the molecules are connected to each reference molecule for which RBFE calculations are quick to converge. The predicted binding free energies are compared to the experimental values in Fig. 8. The precision and accuracy of the NetBFE calculations are summarized in Table 2. In this case, including more reference values does not significantly improve either the precision or the accuracy.

In this work, we have developed a new method to improve the overall statistical precision of binding free energy calculations for multiple molecules against a common biomolecular target. NetBFE (Network Binding Free Energy) implements and extends a previous theoretical work from our group on A-optimal experimental design for measuring pairwise differences, making it practical for realistic binding free energy calculations. We have demonstrated that, at a fixed total computational cost, NetBFE can reduce the statistical errors in BFE calculations by approximately two-fold. The improved precision will increase the confidence in the computational ranking of candidate drug molecules, which is instrumental to prioritizing molecules for synthesis in the design-make-test cycle. In addition, as binding free energy calculations are increasingly used for validating and benchmarking force field models and free energy methods^{4, 11, 7, 50}, the reduced statistical errors afforded by NetBFE will make the benchmark assessments more reliable and thus accelerate the improvements in in silico binding assays.