Statistical Mechanics of Ligand–Receptor Noncovalent Association, Revisited: Binding Site and Standard State Volumes in Modern Alchemical Theories

Methodological and force field effects in the molecular dynamics-based prediction of binding free energies of host-guest systems

As a contribution to the understanding and rationalization of methodological and modeling effects in recent host-guest SAMPL challenges, using an alchemical molecular dynamics technique we have examined the impact of force field parameterization and ionic strength in connection with guest charge neutralization on computed dissociation free energies in two typical SAMPL heavily charged macrocyclic hosts encapsulating small protonated amines with disparate binding affinities. We have shown that the methodological treatment for host neutralization, with explicit ions or with the background neutralizing plasma in the context of alchemical calculations under periodic boundary conditions, has a moderate effect on the calculated affinities. On the other hand, we have shown that seemingly small differences in the force field parameterization in highly symmetric hosts can produce systematic effects on the structural features that can have a significant impact on the predicted binding affinities.

Alchemical Free Energy Estimators and Molecular Dynamics Engines: Accuracy, Precision, and Reproducibility

The binding free energy between a ligand and its target protein is an essential quantity to know at all stages of the drug discovery pipeline. Assessing this value computationally can offer insight into where efforts should be focused in the pursuit of effective therapeutics to treat a myriad of diseases. In this work, we examine the computation of alchemical relative binding free energies with an eye for assessing reproducibility across popular molecular dynamics packages and free energy estimators. The focus of this work is on 54 ligand transformations from a diverse set of protein targets: MCL1, PTP1B, TYK2, CDK2, and thrombin. These targets are studied with three popular molecular dynamics packages: OpenMM, NAMD2, and NAMD3 alpha. Trajectories collected with these packages are used to compare relative binding free energies calculated with thermodynamic integration and free energy perturbation methods. The resulting binding free energies show good agreement between molecular dynamics packages with an average mean unsigned error between them of 0.50 kcal/mol. The correlation between packages is very good, with the lowest Spearman's, Pearson's and Kendall's tau correlation coefficients being 0.92, 0.91, and 0.76, respectively. Agreement between thermodynamic integration and free energy perturbation is shown to be very good when using ensemble averaging.

Relative Binding Free Energy between Chemically Distant Compounds Using a Bidirectional Nonequilibrium Approach

In the context of advanced hit-to-lead drug design based on atomistic molecular dynamics simulations, we propose a dual topology alchemical approach for calculating the relative binding free energy (RBFE) between two chemically distant compounds. The method (termed NE-RBFE) relies on the enhanced sampling of the end-states in bulk and in the bound state via Hamiltonian Replica Exchange, alchemically connected by a series of independent and fast nonequilibrium (NE) simulations. The technique has been implemented in a bidirectional fashion, applying the Crooks theorem to the NE work distributions for RBFE predictions. The dissipation of the NE process, negatively affecting accuracy, has been minimized by introducing a smooth regularization based on shifted electrostatic and Lennard-Jones non bonded potentials. As a challenging testbed, we have applied our method to the calculation of the RBFEs in the recent host–guest SAMPL international contest, featuring a macrocyclic host with guests varying in the net charge, volume, and chemical fingerprints. Closure validation has been successfully verified in cycles involving compounds with disparate Tanimoto coefficients, volume, and net charge. NE-RBFE is specifically tailored for massively parallel facilities and can be used with little or no code modification on most of the popular software packages supporting nonequilibrium alchemical simulations, such as Gromacs, Amber, NAMD, or OpenMM. The proposed methodology bypasses most of the entanglements and limitations of the standard single topology RBFE approach for strictly congeneric series based on free-energy perturbation, such as slowly relaxing cavity water, sampling issues along the alchemical stratification, and the need for highly overlapping molecular fingerprints.

SAMPL9 blind predictions using nonequilibrium alchemical approaches

We present our blind predictions for the Statistical Assessment of the Modeling of Proteins and Ligands (SAMPL), ninth challenge, focusing on the binding of WP6 (carboxy-pillar[6]arene) with ammonium/diammonium cationic guests. Host-guest binding free energies have been calculated using the recently developed virtual double system single box approach, based on the enhanced sampling of the bound and unbound end-states followed by fast switching nonequilibrium alchemical simulations [M. Macchiagodena et al., J. Chem. Theory Comput. 16, 7160 (2020)]. As far as Pearson and Kendall coefficients are concerned, performances were acceptable and, in general, better than those we submitted for calixarenes, cucurbituril-like open cavitand, and beta-cyclodextrines in previous SAMPL host-guest challenges, confirming the reliability of nonequilibrium approaches for absolute binding free energy calculations. In comparison with previous submissions, we found a rather large mean signed error that we attribute to the way the finite charge correction was addressed through the assumption of a neutralizing background plasma.

Seeding the multi-dimensional nonequilibrium pulling for Hamiltonian variation: indirect nonequilibrium free energy simulations at QM levels

The combination of free energy simulations in the alchemical and configurational spaces provides a feasible route to access the thermodynamic profiles under a computationally demanding target Hamiltonian. Normally, due to the significant differences between the computational cost of ab initio quantum mechanics (QM) calculations and those of semi-empirical quantum mechanics (SQM) and molecular mechanics (MM), this indirect method could be used to obtain the QM thermodynamics by combining the SQM or MM results and the SQM-to-QM or MM-to-QM corrections. In our previous work, a multi-dimensional nonequilibrium pulling framework for Hamiltonian variations was introduced based on bidirectional pulling and bidirectional reweighting. The method performs nonequilibrium free energy simulations in the configurational space to obtain the thermodynamic profile along the conformational change pathway under a selected computationally efficient Hamiltonian, and uses the nonequilibrium alchemical method to correct or perturb the thermodynamic profile to that under the target Hamiltonian. The BAR-based method is designed to achieve the best generality and transferability and thus leads to modest (∼20 fold) speedup. In this work, we explore the possibility of further accelerating the nonequilibrium free energy simulation by employing unidirectional pulling and using the selection criterion to obtain the initial configurations used to initiate nonequilibrium trajectories following the idea of adaptive steered molecular dynamics (ASMD). A single initial condition is used to seed the whole multi-dimensional nonequilibrium free energy simulation and the sampling is performed fully in the nonequilibrium ensemble. Introducing very short ps-length equilibrium sampling to grab more initial seeds could also be helpful. The ASMD scheme estimates the free energy difference with the unidirectional exponential average (EXP), but it does not follow exactly the requirements of the EXP estimator. Another deficiency of the seeding simulation is the inherently sequential or serial pulling due to the inter-segment dependency, which triggers some problems in the parallelizability of the simulation. Numerical tests are performed to grasp some insights and guidelines for using this selection-criterion-based ASMD scheme. The presented selection-criterion-based multi-dimensional ASMD scheme follows the same perturbation network of the BAR-based method, and thus could be used in various Hamiltonian-variation cases.

Characterization of the non-covalent interaction between the PF-07321332 inhibitor and the SARS-CoV-2 main protease

We have studied the non-covalent interaction between PF-07321332 and SARS-CoV-2 main protease at the atomic level using a computational approach based on extensive molecular dynamics simulations with explicit solvent. PF-07321332, whose chemical structure has been recently disclosed, is a promising oral antiviral clinical candidate with well-established anti-SARS-CoV-2 activity in vitro. The drug, currently in phase III clinical trials in combination with ritonavir, relies on the electrophilic attack of a nitrile warhead to the catalytic cysteine of the protease. Nonbonded interaction between the inhibitor and the residues of the binding pocket, as well as with water molecules on the protein surface, have been characterized using two different force fields and the two possible protonation states of the main protease catalytic dyad HIS41-CYS145. When the catalytic dyad is in the neutral state, the non-covalent binding is likely to be stronger. Molecular dynamics simulations seems to lend support for an inhibitory mechanism in two steps: a first non-covalent addition with the dyad in neutral form and then the formation of the thiolate-imidazolium ion pair and the ligand relocation for finalising the electrophilic attack.

On the NS-DSSB unidirectional estimates in the SAMPL6 SAMPLing challenge

In the context of the recent SAMPL6 SAMPLing challenge (Rizzi et al. 2020 in J Comput Aided Mol Des 34:601–633) aimed at assessing convergence properties and reproducibility of molecular dynamics binding free energy methodologies, we propose a simple explanation of the severe errors observed in the nonequilibrium switch double-system-single-box (NS-DSSB) approach when using unidirectional estimates. At the same time, we suggest a straightforward and minimal modification of the NS-DSSB protocol for obtaining reliable unidirectional estimates for the process where the ligand is decoupled in the bound state and recoupled in the bulk.

Modeling Kinetics and Thermodynamics of Guest Encapsulation into the [M4L6]12- Supramolecular Organometallic Cage

The encapsulation of molecular guests into supramolecular hosts is a complex molecular recognition process in which the guest displaces the solvent from the host cavity, while the host deforms to let the guest in. An atomistic description of the association would provide valuable insights on the physicochemical properties that guide it. This understanding may be used to design novel host assemblies with improved properties (i.e., affinities) toward a given class of guests. Molecular simulations may be conveniently used to model the association processes. It is thus of interest to establish efficient protocols to trace the encapsulation process and to predict the associated magnitudes ΔG_bind and ΔG_bind^⧧. Here, we report the calculation of the Gibbs energy barrier and Gibbs binding energy by means of explicit solvent molecular simulations for the [Ga₄L₆]^12– metallocage encapsulating a series of cationic molecules. The ΔG_bind^⧧ for encapsulation was estimated by means of umbrella sampling simulations. The steps involved were identified, including ion-pair formation and naphthalene rotation (from L ligands of the metallocage) during the guest’s entrance. The ΔG_bind values were computed using the attach–pull–release method. The results reveal the sensitivity of the estimates on the force field parameters, in particular on atomic charges, showing that higher accuracy is obtained when charges are derived from implicit solvent quantum chemical calculations. Correlation analysis identified some indicators for the binding affinity trends. All computed magnitudes are in very good agreement with experimental observations. This work provides, on one side, a benchmarked way to computationally model a highly charged metallocage encapsulation process. This includes a nonstandard parameterization and charge derivation procedure. On the other hand, it gives specific mechanistic information on the binding processes of [Ga₄L₆]^12– at the molecular level where key motions are depicted. Taken together, the study provides an interesting option for the future design of metal–organic cages.

Binding free energy predictions in host-guest systems using Autodock4. A retrospective analysis on SAMPL6, SAMPL7 and SAMPL8 challenges

We systematically tested the Autodock4 docking program for absolute binding free energy predictions using the host-guest systems from the recent SAMPL6, SAMPL7 and SAMPL8 challenges. We found that Autodock4 behaves surprisingly well, outperforming in many instances expensive molecular dynamics or quantum chemistry techniques, with an extremely favorable benefit-cost ratio. Some interesting features of Autodock4 predictions are revealed, yielding valuable hints on the overall reliability of docking screening campaigns in drug discovery projects.

Binding thermodynamics and interaction patterns of human purine nucleoside phosphorylase-inhibitor complexes from extensive free energy calculations

Human purine nucleoside phosphorylase (hPNP) plays a significant role in the catabolism of deoxyguanosine. The trimeric protein is an important target in the treatment of T-cell cancers and autoimmune disorders. Experimental studies on the inhibition of the hPNP observe that the first ligand bound to one of three subunits effectively inhibits the protein, while the binding of more ligands to the subsequent sites shows negative cooperativities. In this work, we performed extensive end-point and alchemical free energy calculations to determine the binding thermodynamics of the trimeric protein-ligand system. 13 Immucillin inhibitors with experimental results are under calculation. Two widely accepted charge schemes for small molecules including AM1-BCC and RESP are adopted for ligands. The results of RESP are in better agreement with the experimental reference. Further investigations of the interaction networks in the protein-ligand complexes reveal that several residues play significant roles in stabilizing the complex structure. The most commonly observed ones include PHE200, GLU201, MET219, and ASN243. The conformations of the protein in different protein-ligand complexes are observed to be similar. We expect these insights to aid the development of potent drugs targeting hPNP.

SAMPL7 blind predictions using nonequilibrium alchemical approaches

In the context of the SAMPL7 challenge, we computed, employing a non-equilibrium (NE) alchemical technique, the standard binding free energy of two series of host-guest systems, involving as a host the Isaac's TrimerTrip, a Cucurbituril-like open cavitand, and the Gilson's Cyclodextrin derivatives. The adopted NE alchemy combines enhanced sampling molecular dynamics simulations with driven fast out-of-equilibrium alchemical trajectories to recover the free energy via the Jarzynski and Crooks NE theorems. The GAFF2 non-polarizable force field was used for the parametrization. Performances were acceptable and similar in accuracy to those we submitted for Gibb's Deep Cavity Cavitands in the previous SAMPL6 host-guest challenge, confirming the reliability of the computational approach and exposing, in some cases, some important deficiencies of the GAFF2 non-polarizable force field.

Binding Thermodynamics and Interaction Patterns of Inhibitor-Major Urinary Protein-I Binding from Extensive Free-Energy Calculations: Benchmarking AMBER Force Fields

Mouse major urinary protein (MUP) plays a key role in the pheromone communication system. The one-end-closed β-barrel of MUP-I forms a small, deep, and hydrophobic central cavity, which could accommodate structurally diverse ligands. Previous computational studies employed old protein force fields and short simulation times to determine the binding thermodynamics or investigated only a small number of structurally similar ligands, which resulted in sampled regions far from the experimental structure, nonconverged sampling outcomes, and limited understanding of the possible interaction patterns that the cavity could produce. In this work, extensive end-point and alchemical free-energy calculations with advanced protein force fields were performed to determine the binding thermodynamics of a series of MUP-inhibitor systems and investigate the inter- and intramolecular interaction patterns. Three series of inhibitors with a total of 14 ligands were simulated. We independently simulated the MUP-inhibitor complexes under two advanced AMBER force fields. Our benchmark test showed that the advanced AMBER force fields including AMBER19SB and AMBER14SB provided better descriptions of the system, and the backbone root-mean-square deviation (RMSD) was significantly lowered compared with previous computational studies with old protein force fields. Surprisingly, although the latest AMBER force field AMBER19SB provided better descriptions of various observables, it neither improved the binding thermodynamics nor lowered the backbone RMSD compared with the previously proposed and widely used AMBER14SB. The older but widely used AMBER14SB actually achieved better performance in the prediction of binding affinities from the alchemical and end-point free-energy calculations. We further analyzed the protein-ligand interaction networks to identify important residues stabilizing the bound structure. Six residues including PHE38, LEU40, PHE90, ALA103, LEU105, and TYR120 were found to contribute the most significant part of protein-ligand interactions, and 10 residues were found to provide favorable interactions stabilizing the bound state. The two AMBER force fields gave extremely similar interaction networks, and the secondary structures also showed similar behavior. Thus, the intra- and intermolecular interaction networks described with the two AMBER force fields are similar. Therefore, AMBER14SB could still be the default option in free-energy calculations to achieve highly accurate binding thermodynamics and interaction patterns.

Methodological uncertainties in drug-receptor binding free energy predictions based on classical molecular dynamics

Computational approaches are becoming an essential tool in modern drug design and discovery, with fast compound triaging using a combination of machine learning and docking techniques followed by molecular dynamics binding free energies assessment using alchemical techniques. The traditional MD-based alchemical free energy perturbation (FEP) method faces severe sampling issues that may limits its reliability in automated workflows. Here we review the major sources of uncertainty in FEP protocols for drug discovery, showing how the sampling problem can be effectively tackled by switching to nonequilibrium alchemical techniques.

Virtual Double-System Single-Box: A Nonequilibrium Alchemical Technique for Absolute Binding Free Energy Calculations: Application to Ligands of the SARS-CoV-2 Main Protease

In the context of drug-receptor binding affinity calculations using molecular dynamics techniques, we implemented a combination of Hamiltonian replica exchange (HREM) and a novel nonequilibrium alchemical methodology, called virtual double-system single-box, with increased accuracy, precision, and efficiency with respect to the standard nonequilibrium approaches. The method has been applied for the determination of absolute binding free energies of 16 newly designed noncovalent ligands of the main protease (3CL^pro) of SARS-CoV-2. The core structures of 3CL^pro ligands were previously identified using a multimodal structure-based ligand design in combination with docking techniques. The calculated binding free energies for four additional ligands with known activity (either for SARS-CoV or SARS-CoV-2 main protease) are also reported. The nature of binding in the 3CL^pro active site and the involved residues besides the CYS-HYS catalytic dyad have been thoroughly characterized by enhanced sampling simulations of the bound state. We have identified several noncongeneric compounds with predicted low micromolar activity for 3CL^pro inhibition, which may constitute possible lead compounds for the development of antiviral agents in Covid-19 treatment.

SAMPL7 TrimerTrip host-guest binding affinities from extensive alchemical and end-point free energy calculations

The prediction of host-guest binding affinities with computational modelling is still a challenging task. In the 7th statistical assessment of the modeling of proteins and ligands (SAMPL) challenge, a new host named TrimerTrip was synthesized and the thermodynamic parameters of 16 structurally diverse guests binding to the host were characterized. In the TrimerTrip-guest challenge, only structures of the host and the guests are provided, which indicates that the predictions of both the binding poses and the binding affinities are under assessment. In this work, starting from the binding poses obtained from our previous enhanced sampling simulations in the configurational space, we perform extensive alchemical and end-point free energy calculations to calculate the host-guest binding affinities retrospectively. The alchemical predictions with two widely accepted charge schemes (i.e. AM1-BCC and RESP) are in good agreement with the experimental reference, while the end-point estimates perform poorly in reproducing the experimental binding affinities. Aside from the absolute value of the binding affinity, the rank of binding free energies is also crucial in drug design. Surprisingly, the end-point MM/PBSA method seems very powerful in reproducing the experimental rank of binding affinities. Although the length of our simulations is long and the intermediate spacing is dense, the convergence behavior is not very good, which may arise from the flexibility of the host molecule. Enhanced sampling techniques in the configurational space may be required to obtain fully converged sampling. Further, as the length of sampling in alchemical free energy calculations already achieves several hundred ns, performing direct simulations of the binding/unbinding event in the physical space could be more useful and insightful. More details about the binding pathway and mechanism could be obtained in this way. The nonequilibrium method could also be a nice choice if one insists to use the alchemical method, as the intermediate sampling is avoided to some extent.

SAMPL7 TrimerTrip host-guest binding poses and binding affinities from spherical-coordinates-biased simulations

Host-guest binding remains a major challenge in modern computational modelling. The newest 7th statistical assessment of the modeling of proteins and ligands (SAMPL) challenge contains a new series of host-guest systems. The TrimerTrip host binds to 16 structurally diverse guests. Previously, we have successfully employed the spherical coordinates as the collective variables coupled with the enhanced sampling technique metadynamics to enhance the sampling of the binding/unbinding event, search for possible binding poses and calculate the binding affinities in all three host-guest binding cases of the 6th SAMPL challenge. In this work, we report a retrospective study on the TrimerTrip host-guest systems by employing the same protocol to investigate the TrimerTrip host in the SAMPL7 challenge. As no binding pose is provided by the SAMPL7 host, our simulations initiate from randomly selected configurations and are proceeded long enough to obtain converged free energy estimates and search for possible binding poses. The calculated binding affinities are in good agreement with the experimental reference, and the obtained binding poses serve as a nice starting point for end-point or alchemical free energy calculations. Note that as the work is performed after the close of the SAMPL7 challenge, we do not participate in the challenge and the results are not formally submitted to the SAMPL7 challenge.

Upgraded AMBER Force Field for Zinc-Binding Residues and Ligands for Predicting Structural Properties and Binding Affinities in Zinc-Proteins

We developed a novel force field in the context of AMBER parameterization for glutamate and aspartate zinc(II)-binding residues. The interaction between the zinc ion and the coordinating atoms is represented by a spherical nonbonded parameterization. The polarization effect due to the zinc ion has been taken into account by redefining the atomic charges on the residues through accurate quantum mechanical calculations. The new zinc-binding ASP and GLU residues, along with the CYS and HIS zinc-binding residues, parameterized in a recent work [Macchiagodena M.;J. Chem. Inf. Model.2019, 59, 3803-3816], allow users to reliably simulate 96% of the Zn-proteins available in the Protein Data Bank. The upgraded force field for zinc(II)-bound residues has been tested performing molecular dynamics simulations with an explicit solvent and comparing the structural information with experimental data for five different proteins binding zinc(II) with GLU, ASP, HIS, and CYS. We further validated our approach by evaluating the binding free energy of (R)-2-benzyl-3-nitropropanoic acid to carboxypeptidase A using a recently developed nonequilibrium alchemical method. We demonstrated that in this setting it is crucial to take into account polarization effects also on the metal-bound inhibitor.

Identification of potential binders of the main protease 3CLpro of the COVID-19 via structure-based ligand design and molecular modeling

We have applied a computational strategy, using a combination of virtual screening, docking and molecular dynamics techniques, aimed at identifying possible lead compounds for the non-covalent inhibition of the main protease 3CL^pro of the SARS-CoV2 Coronavirus. Based on the X-ray structure (PDB code: 6LU7), ligands were generated using a multimodal structure-based design and then docked to the monomer in the active state. Docking calculations show that ligand-binding is strikingly similar in SARS-CoV and SARS-CoV2 main proteases. The most potent docked ligands are found to share a common binding pattern with aromatic moieties connected by rotatable bonds in a pseudo-linear arrangement.

A remark on the efficiency of the double-system/single-box nonequilibrium approach in the SAMPL6 SAMPLing challenge

The alchemical nonequilibrium switching technique was one of several methods in the top tier of performance in the recent SAMPL6 SAMPLing challenge in both accuracy and efficiency. In this paper, in the context of nonequilibrium alchemical switching, we compare the efficiency of the double-system/single-box (DSSB) approach (used in the SAMPL6 challenges) to the standard single-system/double-box method (SSDB). Exploiting the Crooks theorem in a simple but effective test case, we analytically show that the DSSB approach is almost twice as efficient as SSDB for slow near-equilibrium switching but it gives basically no gain over the conventional SSDB approach when the variance of the work distribution exceeds few [Formula: see text], with the potential of producing artifacts and entanglements if not judiciously implemented.

SAMPL6 host-guest binding affinities and binding poses from spherical-coordinates-biased simulations

Host-guest binding is a challenging problem in computer simulation. The prediction of binding affinities between hosts and guests is an important part of the statistical assessment of the modeling of proteins and ligands (SAMPL) challenges. In this work, the volume-based variant of well-tempered metadynamics is employed to calculate the binding affinities of the host-guest systems in the SAMPL6 challenge. By biasing the spherical coordinates describing the relative position of the host and the guest, the initial-configuration-induced bias vanishes and all possible binding poses are explored. The agreement between the predictions and the experimental results and the observation of new binding poses indicate that the volume-based technique serves as a nice candidate for the calculation of binding free energies and the search of the binding poses.

Upgrading and Validation of the AMBER Force Field for Histidine and Cysteine Zinc(II)-Binding Residues in Sites with Four Protein Ligands

We developed and validated a novel force field in the context of the AMBER parameterization for the simulation of zinc(II)-binding proteins. The proposed force field assumes nonbonded spherical interactions between the central zinc(II) and the coordinating residues. A crucial innovative aspect of our approach is to account for the polarization effects of the cation by redefining the atomic charges of the coordinating residues and an adjustment of Lennard-Jones parameters of Zn-interacting atoms to reproduce mean distance distributions. The optimal transferable parametrization was obtained by performing accurate quantum mechanical calculations on a training set of high-quality protein structures, encompassing the most common folds of zinc(II) sites. The addressed sites contain a zinc(II) ion tetra-coordinated by histidine and cysteine residues and represent about 70% of all physiologically relevant zinc(II) sites in the Protein Data Bank. Molecular dynamics simulations with explicit solvent, carried out on several zinc(II)-binding proteins not included in the training set, show that our model for zinc(II) sites preserves the tetra-coordination of the metal site with remarkable stability, yielding zinc(II)-X mean distances similar to experimental data. Finally, the model was tested by evaluating the zinc(II)-binding affinities, using the alchemical free energy perturbation approach. The calculated dissociation constants correlate satisfactorily with the experimental counterpart demonstrating the validity and transferability of the proposed parameterization for zinc(II)-binding proteins.

Solvation free energies via alchemical simulations: let's get honest about sampling, once more

Free energy perturbation (FEP) approaches with stratification have seen widespread and increasing use in computational studies of biologically relevant molecules. However, when the molecular systems are characterized by a complex conformational free energy landscape, the assessment of convergence remains a concern for many practitioners. The sampling problem in FEP has been authoritatively addressed in a recent perspective paper [D. Mobley, J. Comput.-Aided Mol. Des., 2012, 26, 93], incisively entitled "Let's get honest about sampling". Here, I return to the issue of sampling in the determination of the octanol-water partition coefficient for a synthetic precursor of kinase inhibitors that has been included in the recent extension of the SAMPL6 blind challenge of log P coefficients. I will show that even for this simple compound, whose conformational space is essentially dictated by two sp₃ rotable bonds connecting rigid planar units, canonical sampling using standard techniques can be surprisingly hard to achieve. I will also show how the conformational sampling problem can be effectively bypassed using unidirectional and bidirectional nonequilibrium work methods, reliably recovering the solvation energy with minimal methodological uncertainty.

Understanding PIM-1 kinase inhibitor interactions with free energy simulation

The proviral integration site of the Moloney leukemia virus (PIM) family includes three homologous members. PIM-1 kinase is an important target in effective therapeutic interventions of lymphomas, prostate cancer and leukemia. In the current work, we performed free energy calculations to calculate the binding affinities of several inhibitors targeting this protein. The alchemical method with integration and perturbation-based estimators and the end-point methods were compared. The computational results indicated that the alchemical method can accurately predict the binding affinities, while the end-point methods give relatively unreliable predictions. Decomposing the free energy difference into enthalpic and entropic components with MBAR reweighting enabled us to investigate the detailed thermodynamic parameters with which the entropy-enthalpy compensation in this protein-ligand binding case is identified. We then studied the conformational ensemble, and the important protein-ligand interactions were identified. The current work sheds light on the understanding of the PIM-1-kinase-inhibitor interactions at the atomic level and will be useful in the further development of potential drugs.

BAR-based multi-dimensional nonequilibrium pulling for indirect construction of a QM/MM free energy landscape

Construction of free energy landscapes at the quantum mechanics (QM) level is computationally demanding. As shown in previous studies, by employing an indirect scheme (i.e. constructing a thermodynamic cycle connecting QM states via an alchemical pathway), simulations are converged with much less computational burden. The indirect scheme makes QM/molecular mechanics (MM) free energy simulation orders of magnitude faster than the direct QM/MM schemes. However, the indirect QM/MM simulations were mostly equilibrium sampling based and the nonequilibrium methods were merely exploited in one-dimensional alchemical QM/MM end-state correction at two end states. In this work, we represent a multi-dimensional nonequilibrium pulling scheme for indirect QM/MM free energy simulations, where the whole free energy simulation is performed only with nonequilibrium methods. The collective variable (CV) space we explore is a combination of one alchemical CV and one physically meaningful CV. The current nonequilibrium indirect QM/MM simulation method can be seen as the generalization of equilibrium perturbation based indirect QM/MM methods. The test systems include one backbone dihedral case and one distance case. The two cases are significantly different in size, enabling us to investigate the dependence of the speedup of the indirect scheme on the size of the system. It is shown that the speedup becomes larger when the size of the system becomes larger, which is consistent with the scaling behavior of QM Hamiltonians.

Dissociation Process of a MDM2/p53 Complex Investigated by Parallel Cascade Selection Molecular Dynamics and the Markov State Model

Recently, we efficiently generated dissociation pathways of a protein-ligand complex without applying force bias with parallel cascade selection molecular dynamics (PaCS-MD) and showed that PaCS-MD in combination with the Markov state model (MSM) yielded a binding free energy comparable to experimental values. In this work, we applied the same procedure to a complex of MDM2 protein and the transactivation domain of p53 protein (TAD-p53), the latter of which is known to be very flexible in the unbound state. Using 30 independent MD simulations in PaCS-MD, we successfully generated 25 dissociation pathways of the complex, which showed complete or partial unfolding of the helical region of TAD-p53 during the dissociation process within an average simulation time of 154.8 ± 46.4 ns. The standard binding free energy obtained in combination with one-dimensional-, three-dimensional (3D)- or C_α-MSM was in good agreement with those determined experimentally. Using 3D-MSM based on the center of mass position of TAD-p53 relative to MDM2, the dissociation rate constant was calculated, which was comparable to those measured experimentally. C_α-MSM based on all C_α coordinates of TAD-p53 reproduced the experimentally measured standard binding free energy, and dissociation and association rate constants. We conclude that the combination of PaCS-MD and MSM offers an efficient computational procedure to calculate binding free energies and kinetic rates.