Optimization of Protein–Ligand Electrostatic Interactions Using an Alchemical Free-Energy Method

We present an explicit solvent alchemical free-energy method for optimizing the partial charges of a ligand to maximize the binding affinity with a receptor. This methodology can be applied to known ligand–protein complexes to determine an optimized set of ligand partial atomic changes. Three protein–ligand complexes have been optimized in this work: FXa, P38, and the androgen receptor. The sets of optimized charges can be used to identify design principles for chemical changes to the ligands which improve the binding affinity for all three systems. In this work, beneficial chemical mutations are generated from these principles and the resulting molecules tested using free-energy perturbation calculations. We show that three quarters of our chemical changes are predicted to improve the binding affinity, with an average improvement for the beneficial mutations of approximately 1 kcal/mol. In the cases where experimental data are available, the agreement between prediction and experiment is also good. The results demonstrate that charge optimization in explicit solvent is a useful tool for predicting beneficial chemical changes such as pyridinations, fluorinations, and oxygen to sulfur mutations.

Figure (S1) below shows the three ligands in optimized in this work with each atom named for later reference.

Ligand Charges
Parameters for the ligands in Figure (S1) were generated using Antechamber 4 with AMBER GAFF 2 5 and AM1-BCC 6 . The partial charges for the wild type inhibitors are shown in Table (S1). To inspect the effect of the RMSD on where the optimiser places the charge optimisations were run for 0.01, 0.03 and 0.05 . For each RMSD the optimisation q e was repeated three times the average of the optimised charges across these three replicates for all RMSDs is shown in Table (S2). It can be seen from the low standard deviation for all sets of optimised charges in Table (S2) that the optimisation agrees well on the optimal set of charges across replicates. To visually assert that the charge is being placed in the same places independent of RMSD Figures (S4-S12) show all inhibitors colored by change in charge across all RMSDs. Figures (S4-S12) show that the information about which atoms could be beneficially changed to be more positive or negative is largely invariant as the RMSD is changed.

Convergence Plots
To choose the amount of sampling required to calculate converged the ΔΔG step convergence of the calculation of using SSP for a small perturbation (0.01 ) ΔΔG step e was analyzed. Where is for one step of the optimization algorithm ΔΔG step ΔΔG binding as defined in the methods. To perform this calculation we use the test system FXa and we apply a perturbation of +0.01 to half the atoms and -0.01 to the other half, e e maintaining the simulations net charge. Calculations for the SSP of this ΔΔG step mutant are performed in triplicate scanning across simulation length and presented in Figure (S14).   To compare the speed of SSP and full FEP for computing the gradient of ΔΔG binding w.r.t all charges the convergence of with sampling time needed to be ΔΔG binding investigated when was calculated with full FEP. Figure (S16) presents the ΔΔG binding results of this investigation and shows that for a perturbation of 0.00015 to one e charge of the ligand 1 ns is sufficient to calculated converged . To calculate ΔΔG binding these the protocol is identical to that presented in the FEP Calculations ΔΔG binding section methods, except for the variation in sampling time.

Unconstrained optimization
In this work, a limit is placed on the root mean squared difference (rmsd) between the original and optimized charges on the ligand. Figure (S17, S18, S19) are a presentation of a test of the optimization with no rmsd limit. Figure

Parallelization strategy
The code presented in this work (Ligand charge optimiser) makes an effort to parallelize the work of calculating the dynamics for single step perturbation (SSP), SSP gradient calculations and free energy perturbation (FEP) calculations wherever possible. The decomposition in these three computational tasks is different in each case. Dynamics are decomposed by trajectory length i.e. requesting 50ns over 4 GPUs will compute 12.5 ns per GPU. SSP gradient calculations are decomposed by mutant i.e. requesting the gradient of ddG w.r.t the 16 charges of a ligand over 4 GPUs will compute 4 mutants per GPU. FEP is decomposed by lambda window i.e requesting 16 windows over 4 GPU will compute 4 FEP windows per GPU.