The performance of ensemble-based free energy protocols in computing binding affinities to ROS1 kinase

doi:10.21203/rs.3.rs-1459035/v1

PDF

Article

The performance of ensemble-based free energy protocols in computing binding affinities to ROS1 kinase

https://doi.org/10.21203/rs.3.rs-1459035/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Optimization of binding affinities for compounds to their target protein is a primary objective in drug discovery. Herein we report on a collaborative study that evaluates a set of compounds binding to ROS1 kinase. We use ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent) and TIES (thermodynamic integration with enhanced sampling) protocols to rank the binding free energies. The predicted binding free energies from ESMACS simulations show good correlations with experimental data for subsets of the compounds. Consistent binding free energy differences are generated for TIES and ESMACS. Although an unexplained overestimation exists, we obtain excellent statistical rankings across the set of compounds from the TIES protocol, with a Pearson correlation coefficient of 0.90 between calculated and experimental activities.

ROS1 is a receptor tyrosine kinase closely related to the anaplastic lymphoma kinase (ALK) and leukocyte tyrosine kinase (LTK) based on sequence similarity of their kinase domains. The ROS1 protein is composed of an extracellular domain containing several fibronectin-like repeats and a cytoplasmic kinase domain. Genomic rearrangements involving ROS1 have been detected in a variety of cancers including non-small cell lung cancer (NSCLC), glioblastoma, and colorectal cancer amongst others. Chimeric fusion proteins involving ROS1 fused to the N-terminal domains of different unrelated proteins have been shown to be oncogenic¹. Marketed kinase inhibitor drugs such as Crizotinib and Entrecitinib are active at various kinases including ALK and ROS1 and are used in the treatment of NSCLC. Although the rate of oncogenic ROS1 fusions is generally low, for instance 1–2% in NSCLC, it may be higher in other cancers^2,3. Consequently, the search for ROS1-targeted therapies is the subject of extensive investigation.

The development of a new drug is an expensive and time-consuming process that takes around 10–15 years with an average cost of US $ 1.3–2.9 billion^4,5. Virtual screening approaches have been increasingly applied in the drug discovery process, evidenced by their involvement in more than 70 approved drugs⁶. In most cases, computational approaches have been extensively used in the initial stages including hit to lead and lead optimisation. They can offset the costs associated with bringing novel drugs to market, and accelerate the preclinical drug screening process.

Docking is the most prominently used virtual screening approach. It has been employed as a fast way to predict whether a given compound is able to bind to a specific target, and to estimate the putative bound conformation of the compound. Docking approaches perform well at recapitulating binding modes and assessing the fit of molecules into binding sites. Furthermore, this placement ability allows them to separate active and inactive molecules and provide enriched selections with greater chance of showing bioactivity. However, the docking scores are generally not an accurate indicator for the binding energies and cannot rank the activity of structurally similar molecules. Consequently, more sophisticated approaches are needed that on the one hand capture the dynamic nature of protein ligand binding, whilst also providing a more correct physical assessment of the interaction and binding energies. Molecular dynamics simulations are well suited to study the dynamics of the proteins, refined binding poses, and more importantly, generate more reliable predicted binding affinities.

In recent years, a lot of effort has been invested to develop workflows that simplify and automate the process of free energy calculations⁷, which includes steps to plan, set up and execute simulations, and to analyse the final results. There are a few automated pipelines developed for free energy calculations, including commercial ones like FEP + ⁸ and Molecular Operating Environment (MOE)⁹, and non-commercial ones like Amber free energy workflow (FEW)¹⁰, FESetup¹¹, FEPrepare¹², CHARMM-GUI FEP calculator¹³, to name just a few. The workflows have found promising applications in academia as well as industry. FEP + ⁸, for example, is the free energy calculation suite from Schrödinger Inc. that has made a significant impact in the pharmaceutical industry. To obtain reliable and reproducible free energy estimation, we have developed two ensemble-based protocols in the last few years, named ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent)^7,14 and TIES (thermodynamic integration with enhanced sampling)^7,15. ESMACS is an end-point approach based on the molecular mechanics Poisson-Boltzmann surface area method (MMPBSA)¹⁶, and TIES is an alchemical approach centred on thermodynamic integration (TI). We recently made a systematic uncertainty quantification (UQ) analysis of molecular simulations and binding free energy estimations, and showed that ensembles can be invoked to acquire robust statistical measures of uncertainty, which address both parametric and stochastic components of MD simulations¹⁷. The names of our methods emphasise the central importance of the ensemble-based nature of the protocols employed^7,14,18. Ensemble approaches lead to increased reliability and reproducibility, with tighter control of standard errors revealing the distribution of results that can be obtained^7,19−22. ESMACS and TIES are performed using a binding affinity calculator (BAC)²³, a computational pipeline designed to automate the end-to-end execution of free energy calculations, and to handle ensemble calculations.

The purpose of the present study is to evaluate the ability of ESMACS and TIES to estimate binding affinities of a set of 150 compounds to the ROS1 protein target. Like many protein kinases, the binding site of ROS1 is well structured, and the binding mode for the scaffold of the congeneric series is clearly defined in the crystal structure. Under such circumstances, it is reasonable to expect that good predictions can be achieved, although rotatable bonds in some of the compounds may pose a challenge for conformational sampling and hence convergence of the predictions.

In this section, we first describe the dataset of the compounds and the preparation of the molecular models before explaining the two methodologies used to predict binding affinity values.

Compound set. A set of kinase inhibitor candidates for ROS1 was available from Janssen, which comprises 150 compounds from an internal ROS1 inhibition lead optimization program. The compounds are numbered as per the order in the file provided by Janssen. They are all based on the same chemical scaffold, a tetrahydropyrazolopyrazinone (Fig. 1). The ROS1 activity is provided from a biochemical substrate inhibition assay and a cell proliferation assay. The former is preferred here for comparing with calculated binding free energies. The pIC₅₀ of the compounds ranged from 5.03 to 8.37, a total of 3.34 log units (no activity was detected for five of the compounds at the highest concentration (10 µM) tested; a pIC₅₀ of 5.00 was used hereafter). To capture the experimental variability, we extracted all data for repeat bioactivity measurements of the same compound in the same assay. The average standard deviation in pIC₅₀ for compounds tested at least twice in the ROS1 biochemical assay was 0.12 log units. This is considered an underestimate of the true experimental error which should also include considerations of new compound samples, different assay conditions and laboratories etc.

The compounds differ with functional groups added at different positions of the scaffold (Fig. 1). The compounds can be grouped based on their electrostatic properties (charged and non-charged), or the positions where functional groups are added. The DFG kinase motif is a conserved Asp-Phe-Gly tripeptide at the proximal end of a flexible activation loop. Two distinct conformations, DFG-in and DFG-out, are often observed in the kinase structures, which indicate the catalytic states of the kinases. There are two ROS1 x-ray structures in the protein data bank, both in an active DFG-in conformation. Janssen have obtained four crystallographic structures. All of them have a DFG-out conformation and are complexed with compounds that share the same scaffold (Fig. 1).

Model preparation and molecular simulation. Chain A of one of the structures has been used for the current study. The missing residues in the A-loop are constructed by ModLoop²⁴. All water molecules in the pdb file are retained. There are three phosphorylated tyrosine residues in the x-ray structures; they are phosphorylated as well in the simulations. The compounds were docked with the protein using the OpenEye Posit method²⁵, drawing on information from the bound ligands to improve pose prediction. All the compounds share the same binding pose at the ROS1 binding site, forming one hydrogen bond with the hinge residue (Fig. 2).

Preparation and setup of the simulations were implemented using the BAC (binding affinity calculator)²³, including parametrization of the compounds, solvation of the complexes, electrostatic neutralization of the systems by adding counterions and generation of configurations files for the simulations. Protein parameters were taken from the standard Amber force field for bioorganic systems (ff14SB)²⁶. Ligand parameterizations were produced using the general Amber force field 2 (GAFF2). The partial charges of the compounds were generated using AM1-BCC method which had been suggested to be preferable compared with the RESP method²².

All ligand-protein complexes were solvated in orthorhombic water boxes with a minimum extension from the protein of 14 Å. The protonation states of the residues were assigned using the “reduce” module of AmberTools²⁶. Counterions were added to electrostatically neutralize the systems.

ESMACS. We use the ESMACS¹⁴ protocol for the free energy calculations. We also used TIES¹⁵ for the prediction of relative free energies, described later. The protocols use ensembles of replicas to obtain reproducible binding affinity estimates with robust uncertainty estimates. ESMACS is based on MMPBSA calculations but incorporates a variety of sampling approaches and entropic calculations^27,28. Here conformations of the complex, receptor and compound are all extracted from simulation of the complex, a protocol termed as 1traj²⁹. This is the commonly used approach for the end-point free energy methods. It achieves good convergence because of the cancellation between the noisy terms of the internal energies of the ligand, receptor and complex. Upon binding, however, conformational changes occur for both protein and ligands, which associate with adaptation free energies¹⁴. They are the energy differences between the free state and the bound state. Studies have shown that the inclusion of adaptation energies clearly improves the predictions of binding free energy ranking for some molecular systems^14,27,30, while the situation may be more complicated for others^28,29. In the current study, the binding free energy of multiple compounds are investigated to the same protein receptor. The adaptation energy can be calculated in relative terms using the average energy of the protein²⁷, which is then incorporated into 1traj, a protocol designated as 1traj-ar²⁹.

Ensembles of NAMD³¹ simulations with 25 replicas were performed, with different initial velocities drawn independently from a Maxwell − Boltzmann distribution. A series of equilibration runs, totalling 2 ns, were conducted, with the restraints on heavy atoms gradually removed. Finally, 4 ns production simulations were run for each replica for all ESMACS and TIES simulations. For the ESMACS free energy calculations, the conformations of the protein and the ligands were extracted from the complex simulations.

TIES. We use TIES¹⁵ to predict free energy differences for two groups of compound pairs, using compounds L30 and L2 as the reference structures. The two compounds differ by an isopropyl group which extends into a hydrophobic binding pocket (Fig. 2b). TIES20³² is used to identify the maximum common substructure and to generate the topology and coordinate files for the hybrid compounds. A total of 103 compound pairs have been constructed using the two reference compounds. The compound pairs are chosen so that the atoms in the reference compounds are largely in the common region, with no more than 5 atoms in the alchemical region.

While the first group mainly focuses on the differences at the hydrophobic site, the compounds in the second group share the same isopropyl group at the hydrophobic site but differ at other functional site(s) (Fig. 1). Compound L2 is in both of the groups, and is one of the best binders in the first group. The grouping indeed represents the common lead optimisation steps: optimising one functional position of the lead to find more promising compounds, and then adding/changing functional groups at other positions.

The majority of TIES simulations – 74 out of 103 – are performed with NAMD. The remaining 29 pairs are studied with OpenMM³³. The two MD engines were used to take advantage of the architectures of different supercomputers. We use NAMD on CPU-based ARCHER2 (https://www.archer2.ac.uk/) and SuperMUC-NG (https://doku.lrz.de/display/PUBLIC/SuperMUC-NG), and use OpenMM on the GPU partition of Summit (https://www.olcf.ornl.gov/summit/). Our previous studies have shown that, for the same molecular systems, consistent free energies are obtained not only from different supercomputers using the same MD engine – NAMD in this case³⁰, but also from different MD engines on different supercomputers^21,34 – NAMD³¹ running on CPUs and OpenMM³³ or pmemdGTI³⁵ running on GPUs.

3.1. ESMACS results

A moderate correlation, r = 0.33, is obtained for the entire set of the compounds (Fig. 3). The predicted binding free energies appear to be more negative for the charged compounds. These compounds have no obvious differences with the electrostatically neutral ones in their structural and dynamic behaviours. The inclusion of receptor adaptation energy in the 1traj-ar approach causes the calculated binding free energies to be scattered more widely, with a reduced correlation coefficient (r = 0.22).

The scaffold is aligned well for all the compounds in their initial structures; it has the same pose as shown in Fig. 2. For the compounds with similar experimental binding affinities, the free energies of the charged ones are more favourably predicted than the electrostatically neutral ones (Fig. 3). A good correlation is obtained for the charged compounds, with a correlation coefficient of 0.84. This is largely due to the correct predictions of the ranking for the two compounds with less favourable binding affinities (Fig. 3). The correlation coefficient exhibits no improvements when the charged compounds are excluded, having the same correlation coefficient of 0.33 but a reduced covariance from 1.37 to 1.11 kcal/mol. The weak correlation is due to the wide spread of the predictions for the compounds with less favourable binding affinities. More negative predictions are made for some of these compounds (Fig. 3), which are caused by more negative van der Waals contributions.

Based on the positions and the chemical properties of the modifications made to the scaffold (Fig. 1, Table 1), the compounds are clustered into different groups. The grouping is based on the following categories: (1) ether variations, (2) pyridyl variation, (3) central aryl variation, (4) lactam substituent variation, (5) central aryl and lactam, (6) central aryl and lactam and ethers (Fig. 4). The subgroup within category 6 shows a moderate correlation (r = 0.52), while the subgroup with the combination of 1 and 6 (subgroup 16) manifests a strong correlation (r = 0.66). The combination of the two subgroups, containing about half of the neutral compounds, yields a similar correlation to that for subgroup 16. The other subgroups show weak or very weak correlations/anticorrelations (Fig. 4, Table 1).

3.2. TIES results

The TIES study has been performed for a total of 103 compounds pairs, of which 29 are simulated by OpenMM and 74 by NAMD. Our previous studies have shown that consistent results are generated from different MD engines^21,34 or on different computational platforms³⁰, provided that the same force field is used and the ensemble approach is employed. Here we evaluate the performance of TIES protocols in reproducing the experimental data for the dataset.

Comparison of predictions with experimental data. A very good correlation, with a correlation coefficient of 0.90, is obtained between the calculated free energy differences and the experimental data (Fig. 5a). In the subgroup with compound L30 as the reference, most of the predicted binding free energy differences are negative, indicating that the modifications to the compound are likely to improve the binding potencies (ΔΔG < 0). The calculations correctly predict compound L2 (Fig. 1) as one of the best binders. As L2 already has a very favourable binding free energy (-10.49 kcal/mol), it is not surprising that most of the modifications from it deteriorate the binding affinities (ΔΔG > 0). For all of the 103 compound pairs, 94% of the TIES predictions agree directionally with the experimental observations, that is, the calculated binding free energy differences have the same sign as those from experimental measurements.

The relative binding free energy differences from TIES agree well with those from ESMACS (Fig. 5c), with a correlation coefficient of 0.62. The predicted ΔΔG values for the two subgroups, paired with L30 and L2, are clearly distinct from both of the protocols. The predicted ΔΔG results from ESMACS also have a good correlation with the experimental data, with a correlation coefficient of 0.67 (Fig. 5b).

It should be noted, however, that an apparent overestimation is observed for some of the compound pairs involving L30, making the predicted free energy differences larger than those from experiments. An MSE of -1.36 kcal/mol is obtained for the L30 subgroup, comparing with an MSE of 0.16 kcal/mol for the subgroup with compound L2 as the reference. The apparent overestimations are investigated further in the following sections.

Overestimation from closed cycle mutations. There are 17 compounds which are paired with both L30 and L2 (Fig. 6). This makes it possible to check the performance of the cycle closure method. For the transformations starting from L2 (L2 → L* in Fig. 6a), the average free energy change from TIES calculations agrees well with the experimental data (Fig. 6b). The transformation L30 → L2, however, overestimates the difference by 1.28 kcal/mol. The two-step transformation for all 17 compounds (L30 → L2 → L*) generates an averaged overestimation of 1.17 kcal/mol, similar to that observed in L30 → L2 transformation. Compounds L* and L2 (Fig. 6) have the same/similar alkyl isopropyl substituent at the hydrophobic position (Fig. 1), while L30 does not; transformations L30 → L* therefore share the same change as that in L30 → L2 transformation. There is a similar overestimation in L30 → L* (-0.99 kcal/mol) as that in L30 → L2, which is likely to be rooted in the same alchemical change of growing an alkyl group from L30.

The legs L30 → L2, L2 → L* and L* → L30 generate an average difference of 0.17 kcal/mol within the closed cycle (L30 → L2 → L* → L30, Fig. 6) between the calculations and the experiments. As we have stated previously³⁴, a hysteresis value of 0 is a necessary but not sufficient condition for convergence of predictions. An overall hysteresis near 0 for a closed cycle is not necessarily an indication of convergence for each leg; it may be a result of cancellation from individual legs. It is understandable that extending an alkyl group into a hydrophobic binding pocket (Fig. 2b) makes the binding more favourable, about which TIES and experiment agree with each other directionally. It is likely that the force field and/or the sampling attribute to the differences between the calculations and the experimental measurement, but we cannot pin down the source of the apparent overestimation from the simulations.

Ligand pairs sharing the common change. 50 pairs have been identified which involve a common variation between the two compounds in the pairs, changing from a -CH group to -C-CH-(CH₃)₂ at the hydrophobic site. This is the same change as that between L30 to L2 (Fig. 1). Simulations for the pair get an overestimation of 1.28 kcal/mol for the binding free energy difference: -4.16 kcal/mol from TIES vs -2.88 kcal/mol from experiment. For these 50 pairs, an averaged binding free energy difference of -3.82 kcal/mol is obtained from TIES calculations, while an average of -2.42 kcal/mol from experimental data. A recent study has proposed a combination of machine learning (ML) and alchemical free energy calculation to improve the accuracy of free energy predictions³⁶. The ML method is used to learn the differences between the calculations and the experimental data and derives a correct term which brings the two closer. Such an approach may identify and subsequently correct the discrepancies as observed here. It should be noted that the ML approach used in the study³⁶ only learns the differences between the calculations and experiments, with no attempts to identify the source of the discrepancies. Studies have shown that a large set of typical chemical modifications lead to a nearly normal distribution of binding energy differences³⁷. For individual modifications, however, a wide and non-Gaussian distribution has been observed for experimental binding free energies from a large number of independent measurements³⁸. Such non-normality has also been reported in numerous ensemble simulations when large sample sizes are used^20,22,32. These observations indicate that the assumption of Gaussian statistics may not be justified²⁰. If, on average, these 50 compound pairs have the same discrepancy as that in L30 → L2 transformation, the systematic differences might be collectively “corrected”. If we apply an offset and shift up all of the TIES result by 1.28 kcal/mol, the calculations agree well with the experimental data, with an MSE reducing from − 1.40 kcal/mol to -0.12 kcal/mol (Fig. 7). This offset could be associated with incomplete sampling of the protein conformational movement that is required to accommodate the change from a -CH group to -C-CH-(CH₃)₂ in the region of the DFG loop. DFG loop sampling is known to require longer timescales than studied here³⁹.

Using the TIES and ESMACS protocols, we have computed the binding free energies of a series of ligands interacting with ROS1 kinase. Good correlations are obtained from our ESMACS study for subsets of the compounds. An excellent statistical ranking is generated from the TIES study, with a correlation of 0.90 between the predictions and the experimental measurements. Moreover, TIES and ESMACS agree with one another in the prediction of free energy differences.

An unexplained overestimation is observed in the TIES study for one subset of compound pairs which share a common difference between the pairs. For the transformations to which cycle closure methods can be applied for the error estimation, a cancellation has been observed as the errors are similar for specific legs within the cycle. This results in an overall hysteresis near 0 for a closed cycle, but significant errors remain for individual legs. While continuing improvements in the force fields and the application of statistically robust ensemble methods are essential to generate increased accuracy for the physics-driven alchemical approaches, machine learning methods may help to identify discrepancies between the predictions and experimental data and potentially reduce those differences, providing sufficient quantities of reliable training data are available that incorporate information on the way such data are distributed. Despite the existence of this discrepancy, an excellent correlation is obtained in the free energy differences between the calculations and the experiments, indicating that the TIES approach can be used in prospective drug discovery projects to predict the changes of binding free energies in the lead optimisation process.

Data availability

The structures for the ROS1 protein and the compounds are freely available at http://doi.org/10.23728/b2share.3aa86b110f06482eb4989e819c0270ad

Acknowledgments

We are grateful for funding for the UK MRC Medical Bioinformatics project (grant no. MR/L016311/1), the EPSRC funded UK Consortium on Mesoscale Engineering Sciences (UKCOMES grant no. EP/L00030X/1), the European Commission for EU H2020 CompBioMed2 Centre of Excellence (grant no. 823712) and EU H2020 EXDCI-2 project (grant no. 800957). We acknowledge the Gauss Centre for Supercomputing for providing computing time on the supercomputer SuperMUC-NG (https://doku.lrz.de/display/PUBLIC/SuperMUC-NG) at Leibniz Supercomputing Centre under project COVID-19-SNG1 and the very able assistance of its scientific support staff. We also made use of the ARCHER2 UK National Supercomputing Service (http://www.archer2.ac.uk). Access to ARCHER2 was provided through the UKCOMES grant. Some of this work was performed on Summit supercomputers (https://www.olcf.ornl.gov/summit/) at the Oak Ridge National Laboratory, thanks to a 2021 DOE INCITE award “COMPBIO”.

Author contributions

S.W.: Conceptualization, Data curation, Formal analysis, Methodology, Writing—original draft, review & editing. A.P.B.: Data curation. D.W.W.: Conceptualization, Data curation. A.D.W.: Data curation. G.T.: Conceptualization, Data curation, Writing—review & editing. H.v.V.: Conceptualization, Writing—review & editing. P.V.C.: Funding acquisition, Conceptualization, Methodology, Supervision, Writing—original draft, review & editing.

Competing interests

The authors declare no competing interests.

Drilon, A. et al. ROS1-dependent cancers - biology, diagnostics and therapeutics. Nat Rev Clin Oncol 18, 35–55, doi:10.1038/s41571-020-0408-9 (2021).
Gu, T. L. et al. Survey of tyrosine kinase signaling reveals ROS kinase fusions in human cholangiocarcinoma. PLoS One 6, e15640, doi:10.1371/journal.pone.0015640 (2011).
Wiesner, T. et al. Kinase fusions are frequent in Spitz tumours and spitzoid melanomas. Nat Commun 5, 3116, doi:10.1038/ncomms4116 (2014).
DiMasi, J. A., Grabowski, H. G. & Hansen, R. W. Innovation in the pharmaceutical industry: New estimates of R&D costs. J Health Econ 47, 20–33, doi:10.1016/j.jhealeco.2016.01.012 (2016).
Wouters, O. J., McKee, M. & Luyten, J. Estimated Research and Development Investment Needed to Bring a New Medicine to Market, 2009–2018. JAMA 323, 844–853, doi:10.1001/jama.2020.1166 (2020).
Sabe, V. T. et al. Current trends in computer aided drug design and a highlight of drugs discovered via computational techniques: A review. Eur J Med Chem 224, 113705, doi:10.1016/j.ejmech.2021.113705 (2021).
Wan, S., Bhati, A. P., Zasada, S. J. & Coveney, P. V. Rapid, accurate, precise and reproducible ligand-protein binding free energy prediction. Interface Focus 10, 20200007, doi:10.1098/rsfs.2020.0007 (2020).
Wang, L. 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, doi:10.1021/ja512751q (2015).
Molecular Operating Environment (MOE), 2019.01. Chemical Computing Group ULC, 1010 Sherbooke St. West, Suite #910, Montreal, QC, Canada, H3A 2R7, 2021.
Homeyer, N. & Gohlke, H. FEW: a workflow tool for free energy calculations of ligand binding. J Comput Chem 34, 965–973, doi:10.1002/jcc.23218 (2013).
Loeffler, H. H., Michel, J. & Woods, C. FESetup: Automating Setup for Alchemical Free Energy Simulations. J Chem Inf Model 55, 2485–2490, doi:10.1021/acs.jcim.5b00368 (2015).
Zavitsanou, S. et al. FEPrepare: A Web-Based Tool for Automating the Setup of Relative Binding Free Energy Calculations. J Chem Inf Model 61, 4131–4138, doi:10.1021/acs.jcim.1c00215 (2021).
Kim, S. et al. CHARMM-GUI Free Energy Calculator for Absolute and Relative Ligand Solvation and Binding Free Energy Simulations. J Chem Theory Comput 16, 7207–7218, doi:10.1021/acs.jctc.0c00884 (2020).
Wan, S., Knapp, B., Wright, D. W., Deane, C. M. & Coveney, P. V. Rapid, precise, and reproducible prediction of peptide-MHC binding affinities from molecular dynamics that correlate well with experiment. J Chem Theory Comput 11, 3346–3356, doi:10.1021/acs.jctc.5b00179 (2015).
Bhati, A. P., Wan, S., Wright, D. W. & Coveney, P. V. Rapid, Accurate, Precise, and Reliable Relative Free Energy Prediction Using Ensemble Based Thermodynamic Integration. J Chem Theory Comput 13, 210–222, doi:10.1021/acs.jctc.6b00979 (2017).
Kollman, P. A. et al. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc Chem Res 33, 889–897 (2000).
Vassaux, M., Wan, S., Edeling, W. & Coveney, P. V. Ensembles Are Required to Handle Aleatoric and Parametric Uncertainty in Molecular Dynamics Simulation. J Chem Theory Comput 17, 5187–5197, doi:10.1021/acs.jctc.1c00526 (2021).
Coveney, P. V. & Wan, S. On the calculation of equilibrium thermodynamic properties from molecular dynamics. Phys Chem Chem Phys 18, 30236–30240, doi:10.1039/c6cp02349e (2016).
Knapp, B., Ospina, L. & Deane, C. M. Avoiding False Positive Conclusions in Molecular Simulation: The Importance of Replicas. J Chem Theory Comput 14, 6127–6138, doi:10.1021/acs.jctc.8b00391 (2018).
Wan, S., Sinclair, R. C. & Coveney, P. V. Uncertainty quantification in classical molecular dynamics. Philos Trans R Soc A 379, 20200082, doi:10.1098/rsta.2020.0082 (2021).
Wade, A., Bhati, A. P., Wan, S. & Coveney, P. V. Alchemical Free Energy Estimators and Molecular Dynamics Engines: Accuracy, Precision and Reproducibility. ChemRxiv, DOI: 10.26434/chemrxiv-22021-nqp26438r, doi:10.26434/chemrxiv-2021-nqp8r (2021).
Bhati, A. P. & Coveney, P. V. Large Scale Study of Ligand-Protein Relative Binding Free Energy Calculations: Actionable Predictions from Statistically Robust Protocols. ChemRxiv, DOI: 10.26434/chemrxiv-22021-zdzng, doi:10.26434/chemrxiv-2021-zdzng (2021).
Sadiq, S. K. et al. Automated molecular simulation based binding affinity calculator for ligand-bound HIV-1 proteases. J Chem Inf Model 48, 1909–1919, doi:10.1021/ci8000937 (2008).
Fiser, A. & Sali, A. ModLoop: automated modeling of loops in protein structures. Bioinformatics 19, 2500–2501, doi:10.1093/bioinformatics/btg362 (2003).
Kelley, B. P., Brown, S. P., Warren, G. L. & Muchmore, S. W. POSIT: Flexible Shape-Guided Docking For Pose Prediction. J Chem Inf Model 55, 1771–1780, doi:10.1021/acs.jcim.5b00142 (2015).
Amber 2021. University of California, San Francisco.
Wan, S. et al. Evaluation and Characterization of Trk Kinase Inhibitors for the Treatment of Pain: Reliable Binding Affinity Predictions from Theory and Computation. J Chem Inf Model 57, 897–909, doi:10.1021/acs.jcim.6b00780 (2017).
Wright, D. W. et al. Application of the ESMACS Binding Free Energy Protocol to a Multi-Binding Site Lactate Dehydogenase A Ligand Dataset. Adv Theory Simul 3, 1900194, doi:10.1002/adts.201900194 (2019).
Wright, D. W. et al. Application of ESMACS binding free energy protocols to diverse datasets: Bromodomain-containing protein 4. Sci Rep 9, 6017, doi:10.1038/s41598-019-41758-1 (2019).
Wan, S. et al. Rapid and Reliable Binding Affinity Prediction of Bromodomain Inhibitors: A Computational Study. J Chem Theory Comput 13, 784–795, doi:10.1021/acs.jctc.6b00794 (2017).
Phillips, J. C. et al. Scalable molecular dynamics with NAMD. J Comput Chem 26, 1781–1802, doi:10.1002/jcc.20289 (2005).
Bieniek, M. K., Bhati, A. P., Wan, S. & Coveney, P. V. TIES 20: Relative Binding Free Energy with a Flexible Superimposition Algorithm and Partial Ring Morphing. J Chem Theory Comput 17, 1250–1265, doi:10.1021/acs.jctc.0c01179 (2021).
Eastman, P. et al. OpenMM 7: Rapid development of high performance algorithms for molecular dynamics. PLoS Comput Biol 13, e1005659, doi:10.1371/journal.pcbi.1005659 (2017).
Bhati, A. P., Wan, S., Hu, Y., Sherborne, B. & Coveney, P. V. Uncertainty Quantification in Alchemical Free Energy Methods. J Chem Theory Comput 14, 2867–2880, doi:10.1021/acs.jctc.7b01143 (2018).
Lee, T. S., Hu, Y., Sherborne, B., Guo, Z. & York, D. M. Toward Fast and Accurate Binding Affinity Prediction with pmemdGTI: An Efficient Implementation of GPU-Accelerated Thermodynamic Integration. J Chem Theory Comput 13, 3077–3084, doi:10.1021/acs.jctc.7b00102 (2017).
Scheen, J. et al. Hybrid Alchemical Free Energy/Machine-Learning Methodology for the Computation of Hydration Free Energies. J Chem Inf Model 60, 5331–5339, doi:10.1021/acs.jcim.0c00600 (2020).
Hajduk, P. J. & Sauer, D. R. Statistical analysis of the effects of common chemical substituents on ligand potency. J Med Chem 51, 553–564, doi:10.1021/jm070838y (2008).
Wan, S. et al. Evaluation and Characterization of Isoxazole Amides as SMYD3 Inhibitors. ChemRxiv, doi:10.26434/chemrxiv-2022-hrxqs (2022).
Meng, Y., Lin, Y. L. & Roux, B. Computational study of the "DFG-flip" conformational transition in c-Abl and c-Src tyrosine kinases. J Phys Chem B 119, 1443–1456, doi:10.1021/jp511792a (2015).

Table 1

Correlations of calculated and experimental binding free energies for the subgroups of the compounds.
Subgroup*	0	16	2	3	36	456	56	6	6 + 16
No. of comp	1	36	7	2	23	26	14	31	67
correlation	-	0.66	0.28	-1	0.05	-0.39	-0.37	0.52	0.68

*The subgroups are defined using single categories (1 – 6) or combinations of these categories. 16, for example, is a subgroup of compounds with modifications of (1) ether variations and (6) central aryl and lactam and ethers. “6+16” is a subset of compounds from subgroups 6 and 16. L30 is labelled as the only member in subgroup 0. Only electrostatically neutral compounds are included.

No competing interests reported.

PDF

Editorial decision: Major revision
05 Apr, 2022
Reviews received at journal
29 Mar, 2022
Reviewers agreed at journal
17 Mar, 2022
Reviewers invited by journal
17 Mar, 2022
Editor assigned by journal
17 Mar, 2022
Editor invited by journal
17 Mar, 2022
Submission checks completed at journal
17 Mar, 2022
First submitted to journal
16 Mar, 2022

You are reading this latest preprint version

The performance of ensemble-based free energy protocols in computing binding affinities to ROS1 kinase

Status:

Version 1

Abstract

Figures

1. Introduction

2. Methods And Simulations

3. Results

4. Conclusion

Declarations

References

Tables

Competing Interests

Status:

Version 1