The importance of protonation and tautomerization in relative binding affinity prediction: a comparison of AMBER TI and Schrödinger FEP

Structure preparation of FXa data set

The structures of ligand and receptor were carefully prepared based on the high resolution Factor Xa–L51a crystal complex structure (PDB code: 2RA0, resolution: 2.3Å) with no missing residues. The structure of Factor Xa was first separated from the crystal complex with all crystallographic water molecules. Back mutation was modeled for L88V mutation in the L chain. N-methyl-capping group (NME) was added to the C-termini of both chain A and L, while Acetyl-capping group (ACE) was only added to the N-termini of the chain L which is in the exterior of the protein. The missing but structurally important Ca²⁺ and Na⁺ ions were placed and aligned with PDB structure 2W26. The structure of L51a in the crystal complex structure was used as the starting point to build models for ligand L51b–L51k. For details of the system preparation, please refer to Homeyer and Gohlke’s work in the Ref. [14]. The prepared structures were used as input structures for FEW and Schrödinger FEP+ workflows after checking the correctness by visual inspection on the assignment of residue protonation states, rotamers, disulfide bond connections, ions, termini, ligand atom types and starting conformations. The apparent pKa prediction (pKa^app) from the ACD Labs/pKa DB (v.11.02, Advanced Chemistry Development, Inc., Toronto, Ontario, Canada) is based on experimental data and is widely accepted as providing quality pKa predictions. We estimated the protonated state structures of the inhibitors using the ACD Labs/pKa DB calculation algorithm [16, 17] (see Table 1), leading to significant changes from the neutral state of the inhibitors exclusively used in the original FEW paper [14]. The preferred, protonated input structures of receptor and inhibitors were used in comparisons of the FEP+ and AMBER TI predictions later.

Table 1.

Structures, apparent pKa, experimental K_i and converted binding affinities of the Factor X a inhibitors

Ligand	R1 (Charged)	pKaapp (charged)	Kia	ΔGb
L51a		7.84	15.9	−10.71
L51b		8.41	1.4	−12.15
L51btc		7.55	1.4	−12.15
L51c		8.60	8.1	−11.11
L51d		7.66	2255.0	−7.75
L51ed		7.38	147.0	−9.38
L51f		7.66	3.0	−11.70
L51g		7.42	4.4	−11.47
L51hd		6.91	4.1	−11.51
L51i		8.92	690.0	−8.46
L51k		8.62	4.9	−11.41

Open in a new tab

^a

Experimentally measured K_i (unit: nM) are taken from Ref. [18]

^b

Free energies (unit: kcal/mol) are converted based on experimental K_i values by ΔG = −RT lnK_i with T = 300 K

^c

L51bt is the tautomer structure of L51b

^d

Protonated structure was used as the corrected ligand structure of L51e, L51h since their pKa^app are close to 7.4

AMBER FEW TI calculations

The AMBER FEW thermodynamic integration approach was used to calculate relative binding free energies using the AMBER14 suite [19]. Nine ligand transformations based on the structural similarity score in the Ref. [14] were preserved in this work. Within the workflow, atom types of the inhibitors were automatically assigned from General Amber Force Field (GAFF) [20], atomic partial charges were prepared using the AM1-BCC approach [21, 22], and the relative binding free energies (ΔΔG_bind) between two structurally similar ligands was determined by the dual topology, soft-core approach [23]. For each transformation, the alchemical space was divided into nine λ-values from 0.1 to 0.9 with Δλ = 0.1. Simulations at the endpoints (λ = 0 and λ = 1) were not conducted, as recommended for the soft-core approach in AMBER. A default of 5 ns simulation length was run in each of the nine λ windows. Both of the ff99SB and ff12SB force fields were employed to calculate Δ3ΔG_bind of the original data sets, and the ff12SB were applied as the default force field. The convergence was repeatedly measured every 250 ps, and the simulation at the given λ step will be terminated once it is below the defined threshold.

The differences in the free energies for the complex and ligand transformation were calculated according to Eq. (1), respectively. Binding free energy was computed by numerical integration employing the calculation approaches with and without linear extrapolation of dV/dλ to the physical end states at λ = 0 and 1.

Since the same receptor was used, by taking the difference of these values, Eq. (2) gives the relative binding free energy. (Details of the TI setup protocol and determination of the ligand pairs for the transformation simulations are described in Ref. [14]).

Schrödinger FEP calculations

The relative binding free energies were estimated using the Desmond FEP+ module in Schrödinger Maestro Suite 2015 [24]. The OPLS2005 force field was used to describe the protein and ligand. The setup of the FEP [25] calculation was automated by the FEP+ module, in which the CM1A-BCC method was used to prepare ligand atomic partial charges [26]. Unlike the AMBER TI workflow, the Desmond FEP+ has been implemented on a graphics processing unit (GPU) platform. A FEP/REST (free energy perturbation/replica exchange with solute tempering) algorithm [27, 28] was employed to accelerate conformational sampling by effectively locally heating the binding region [28]. The ligand transformations were determined by using the FEP Mapper module based on the similarity scores. The direct transformations (no intermediate ligands involved) which contain the same ligand pairs as the TI calculations were selected for the final FEP calculations. Cycle closure calculations [29] were used for error analysis. Details of the method were described in Ref. [2].

Re-examination of AMBER force fields, convergence methods and linear extrapolation

In the showcase example of Homeyer and Gohlke’s work [14], the AMBER ff99SB force field was used to conduct the TI calculations for the ligands in their neutral states using convergence method 1 (details of the convergence methods are described in below). By taking the FEW TI workflow approach, we repeated the TI calculations with the same settings and additionally carried out calculations with the newer ff12SB force field. Similar TI ΔΔG predictions were generated as the Ref. [14] using the ff99SB force field with a correlation R² = 0.88 to the result with linear extrapolation at the endpoints (see data in Table S1, Fig. S1). Compared to the AMBER ff99SB, the ff12SB force field markedly improved the correlation of the relative binding free energy with the experimental result regardless of which convergence method was used, shown in Table 2. The ff12SB force field was thereafter used as the default.

Two approaches were employed to check the convergence of the dV/dλ values as implemented in FEW, where method 1 evaluated the difference in the standard error of the mean dV/dλ values [23] with error limit 0.01 kcal/mol in two successive checking steps, and method 2 examined the deviation of the true mean of dV/dλ from the sample mean by using the Student’s distribution and a confidence limit of 95 % with error limit 0.2 kcal/mol. We found method 2 converged slightly more slowly than method 1. However, the differences in the uncertainties of the calculated ΔΔG_bind values are trivial, and there is almost no difference in the overall correlations to the experimental result (details of the result are shown Table S1).

Consistent with the finding in the Ref. [14], the ΔΔG predictions without extrapolation at the end states using both convergence methods were markedly improved over the linear extrapolation result, which is commonly understood by the nonlinear dV/dλ curves.

Improvement of AMBER TI predictions with ligand protonation state and tautomerization correction

In drug discovery, protonation states and tautomerization are easily overlooked. By using the commonly applied apparent pKa prediction tool ACD Labs/pKa DB, of which prediction algorithm [16, 17] is based on experimental data and is widely accepted as providing quality pKa predictions, we determined the protonation states of each ligand at pH = 7.4 (data are shown in Table 1). TI predictions of the relative binding affinities of charged ligand transformations are listed in Table S2. Generally, a thermodynamically stable tautomer is favored at equilibrium. Experimentally, ligand L51b is about 11 fold more potent than L51a, resulting in a −1.45 kcal/mol experimental ΔΔG, while the calculation suggested a +1.52 kcal/mol ΔΔG. Quantum Mechanics calculations indicated the tautomer state of ligand L51b is preferred. By considering the alternative tautomer of ligand L51b (referred as ligand L51bt), we found that the predicted relative binding affinity of L51bt to L51a was significantly improved and ΔΔG was now −1.17 kcal/mol. The rank was false predicted by using L51b structure, where 2.23 kcal/mol (3.73 kcal/mol by using ff99SB) was found by using neutral state of the inhibitors, and 1.52 kcal/mol in charged state; by using the L51bt structure, the ΔΔG_bind reduced to 0.16 kcal/mol in neutral state, and moved to −1.17 kcal/mol in charged state which captured the correct rank in a tolerated error.

The effect on overall correlations of AMBER TI predicted relative binding affinities are shown in Table 2, demonstrating the importance of incorporating ligand state at pH = 7.4 and tautomerization correction. After correcting the ligand structures, the correlation R² of AMBER TI prediction was improved remarkably from 0.41 to 0.63 (up to 0.73 without linear extrapolation), and the root mean squared error (RMSE) was reduced to less than 1.0 kcal/-mol within the experimental error (shown in Fig. 1).

Schrödinger FEP prediction

By taking the Schrödinger FEP+ workflow using Desmond, we conducted the FEP calculations for the ligand transformations determined by the FEP Mapper tool. The uncorrected FEP-predicted relative binding affinities and cycle closure corrected energetics for the ligands in both neutral state and charged state are shown in Table S3. By fitting to the experimental results, a correlation coefficient R² of 0.38 was found for the neutral ligand pair transformations and 0.49 for the charged ligand pairs. Interestingly, the cycle closure correction contributed no improvement to the correlation.

Comparison of Schrödinger FEP and AMBER TI workflows

By taking the same input structures, the calculation results by using Schrödinger FEP (corrected predictions) were highly correlated with the AMBER TI predictions (with linear extrapolation). They are almost equivalent with a correlation R² = 0.80, RMSE = 0.64 kcal/mol at neutral state, and R² = 0.96, RMSE = 0.30 kcal/mol when the charges and protonation states are corrected for all the ligands (the correlation of Amber TI and Schrödinger FEP are shown in Fig. 2).

The Schrödinger FEP and AMBER TI workflows are then comparable except for the speed: for the AMBER TI workflow, it takes approximately 1 week to perform one transformation with TI calculation on a state-of-the-art computer cluster using 64 CPU cores (16 cores/node) per λ window, but it only takes a day or less to complete one Schrödinger FEP calculation using 4 GPU cores per transformation. GPU supported AMBER TI module is in active development and is expected to be available in the AMBER16 release [30].