Resolving coupled pH titrations using alchemical free energy calculations.

In a protein, nearby titratable sites can be coupled: the (de)protonation of one may affect the other. The degree of this interaction depends on several factors and can influence the measured p K a . Here, we derive a formalism based on double free energy differences ( Δ Δ G ) for quantifying the individual site p K a values of coupled residues. As Δ Δ G values can be obtained by means of alchemical free energy calculations, the presented approach allows for a convenient estimation of coupled residue p K a s in practice. We demonstrate that our approach and a previously proposed microscopic p K a formalism, can be combined with alchemical free energy calculations to resolve pH-dependent protein p K a values. Toy models and both, regular and constant-pH molecular dynamics simulations, alongside experimental data, are used to validate this approach. Our results highlight the insights gleaned when coupling and microstate probabilities are analyzed and suggest extensions to more complex enzymatic contexts. Furthermore, we find that naïvely computed p K a values that ignore coupling, can be significantly improved when coupling is accounted for, in some cases reducing the error by half. In short, alchemical free energy methods can resolve the p K a values of both uncoupled and coupled residues.

Current State of Open Source Force Fields in Protein-Ligand Binding Affinity Predictions.

In drug discovery, the in silico prediction of binding affinity is one of the major means to prioritize compounds for synthesis. Alchemical relative binding free energy (RBFE) calculations based on molecular dynamics (MD) simulations are nowadays a popular approach for the accurate affinity ranking of compounds. MD simulations rely on empirical force field parameters, which strongly influence the accuracy of the predicted affinities. Here, we evaluate the ability of six different small-molecule force fields to predict experimental protein-ligand binding affinities in RBFE calculations on a set of 598 ligands and 22 protein targets. The public force fields OpenFF Parsley and Sage, GAFF, and CGenFF show comparable accuracy, while OPLS3e is significantly more accurate. However, a consensus approach using Sage, GAFF, and CGenFF leads to accuracy comparable to OPLS3e. While Parsley and Sage are performing comparably based on aggregated statistics across the whole dataset, there are differences in terms of outliers. Analysis of the force field reveals that improved parameters lead to significant improvement in the accuracy of affinity predictions on subsets of the dataset involving those parameters. Lower accuracy can not only be attributed to the force field parameters but is also dependent on input preparation and sampling convergence of the calculations. Especially large perturbations and nonconverged simulations lead to less accurate predictions. The input structures, Gromacs force field files, as well as the analysis Python notebooks are available on GitHub.

Pre-Exascale Computing of Protein-Ligand Binding Free Energies with Open Source Software for Drug Design.

Nowadays, drug design projects benefit from highly accurate protein-ligand binding free energy predictions based on molecular dynamics simulations. While such calculations have been computationally expensive in the past, we now demonstrate that workflows built on open source software packages can efficiently leverage pre-exascale computing resources to screen hundreds of compounds in a matter of days. We report our results of free energy calculations on a large set of pharmaceutically relevant targets assembled to reflect industrial drug discovery projects.

GROMACS in the Cloud: A Global Supercomputer to Speed Up Alchemical Drug Design.

We assess costs and efficiency of state-of-the-art high-performance cloud computing and compare the results to traditional on-premises compute clusters. Our use case is atomistic simulations carried out with the GROMACS molecular dynamics (MD) toolkit with a particular focus on alchemical protein-ligand binding free energy calculations. We set up a compute cluster in the Amazon Web Services (AWS) cloud that incorporates various different instances with Intel, AMD, and ARM CPUs, some with GPU acceleration. Using representative biomolecular simulation systems, we benchmark how GROMACS performs on individual instances and across multiple instances. Thereby we assess which instances deliver the highest performance and which are the most cost-efficient ones for our use case. We find that, in terms of total costs, including hardware, personnel, room, energy, and cooling, producing MD trajectories in the cloud can be about as cost-efficient as an on-premises cluster given that optimal cloud instances are chosen. Further, we find that high-throughput ligand-screening can be accelerated dramatically by using global cloud resources. For a ligand screening study consisting of 19 872 independent simulations or ∼200 µs of combined simulation trajectory, we made use of diverse hardware available in the cloud at the time of the study. The computations scaled-up to reach peak performance using more than 4 000 instances, 140 000 cores, and 3 000 GPUs simultaneously. Our simulation ensemble finished in about 2 days in the cloud, while weeks would be required to complete the task on a typical on-premises cluster consisting of several hundred nodes.

Non-equilibrium approach for binding free energies in cyclodextrins in SAMPL7: force fields and software.

In the current work we report on our participation in the SAMPL7 challenge calculating absolute free energies of the host-guest systems, where 2 guest molecules were probed against 9 hosts-cyclodextrin and its derivatives. Our submission was based on the non-equilibrium free energy calculation protocol utilizing an averaged consensus result from two force fields (GAFF and CGenFF). The submitted prediction achieved accuracy of [Formula: see text] in terms of the unsigned error averaged over the whole dataset. Subsequently, we further report on the underlying reasons for discrepancies between our calculations and another submission to the SAMPL7 challenge which employed a similar methodology, but disparate ligand and water force fields. As a result we have uncovered a number of issues in the dihedral parameter definition of the GAFF 2 force field. In addition, we identified particular cases in the molecular topologies where different software packages had a different interpretation of the same force field. This latter observation might be of particular relevance for systematic comparisons of molecular simulation software packages. The aforementioned factors have an influence on the final free energy estimates and need to be considered when performing alchemical calculations.

Non-active site mutants of HIV-1 protease influence resistance and sensitisation towards protease inhibitors.

BACKGROUND: HIV-1 can develop resistance to antiretroviral drugs, mainly through mutations within the target regions of the drugs. In HIV-1 protease, a majority of resistance-associated mutations that develop in response to therapy with protease inhibitors are found in the protease's active site that serves also as a binding pocket for the protease inhibitors, thus directly impacting the protease-inhibitor interactions. Some resistance-associated mutations, however, are found in more distant regions, and the exact mechanisms how these mutations affect protease-inhibitor interactions are unclear. Furthermore, some of these mutations, e.g. N88S and L76V, do not only induce resistance to the currently administered drugs, but contrarily induce sensitivity towards other drugs. In this study, mutations N88S and L76V, along with three other resistance-associated mutations, M46I, I50L, and I84V, are analysed by means of molecular dynamics simulations to investigate their role in complexes of the protease with different inhibitors and in different background sequence contexts. RESULTS: Using these simulations for alchemical calculations to estimate the effects of mutations M46I, I50L, I84V, N88S, and L76V on binding free energies shows they are in general in line with the mutations' effect on [Formula: see text] values. For the primary mutation L76V, however, the presence of a background mutation M46I in our analysis influences whether the unfavourable effect of L76V on inhibitor binding is sufficient to outweigh the accompanying reduction in catalytic activity of the protease. Finally, we show that L76V and N88S changes the hydrogen bond stability of these residues with residues D30/K45 and D30/T31/T74, respectively. CONCLUSIONS: We demonstrate that estimating the effect of both binding pocket and distant mutations on inhibitor binding free energy using alchemical calculations can reproduce their effect on the experimentally measured [Formula: see text] values. We show that distant site mutations L76V and N88S affect the hydrogen bond network in the protease's active site, which offers an explanation for the indirect effect of these mutations on inhibitor binding. This work thus provides valuable insights on interplay between primary and background mutations and mechanisms how they affect inhibitor binding.

Polycation-Anionic Lipid Membrane Interactions.

Natural or synthetic polycations are used as biocides or as drug/gene carriers. Understanding the interactions between these macromolecules and cell membranes at the molecular level is therefore of great importance for the design of effective polymer biocides or biocompatible polycation-based delivery systems. Until now, details of the processes at the interface between polycations and biological systems have not been fully recognized. In this study, we consider the effect of strong polycations with quaternary ammonium groups on the properties of anionic lipid membranes that we use as a model system for protein-free cell membranes. For this purpose, we employed experimental measurements and atomic-scale molecular dynamics (MD) simulations. MD simulations reveal that the polycations are strongly hydrated in the aqueous phase and do not lose the water shell after adsorption at the bilayer surface. As a result of strong hydration, the polymer chains reside at the phospholipid headgroup and do not penetrate to the acyl chain region. The polycation adsorption involves the formation of anionic lipid-rich domains, and the density of anionic lipids in these domains depends on the length of the polycation chain. We observed the accumulation of anionic lipids only in the leaflet interacting with the polymer, which leads to the formation of compositionally asymmetric domains. Asymmetric adsorption of the polycation on only one leaflet of the anionic membrane strongly affects the membrane properties in the polycation-membrane contact areas: (i) anionic lipid accumulates in the region near the adsorbed polymer, (ii) acyl chain ordering and lipid packing are reduced, which results in a decrease in the thickness of the bilayer, and (iii) polycation-anionic membrane interactions are strongly influenced by the presence and concentration of salt. Our results provide an atomic-scale description of the interactions of polycations with anionic lipid bilayers and are fully supported by the experimental data. The outcomes are important for understanding the correlation of the structure of polycations with their activity on biomembranes.

The SAMPL6 SAMPLing challenge: assessing the reliability and efficiency of binding free energy calculations.

Approaches for computing small molecule binding free energies based on molecular simulations are now regularly being employed by academic and industry practitioners to study receptor-ligand systems and prioritize the synthesis of small molecules for ligand design. Given the variety of methods and implementations available, it is natural to ask how the convergence rates and final predictions of these methods compare. In this study, we describe the concept and results for the SAMPL6 SAMPLing challenge, the first challenge from the SAMPL series focusing on the assessment of convergence properties and reproducibility of binding free energy methodologies. We provided parameter files, partial charges, and multiple initial geometries for two octa-acid (OA) and one cucurbit[8]uril (CB8) host-guest systems. Participants submitted binding free energy predictions as a function of the number of force and energy evaluations for seven different alchemical and physical-pathway (i.e., potential of mean force and weighted ensemble of trajectories) methodologies implemented with the GROMACS, AMBER, NAMD, or OpenMM simulation engines. To rank the methods, we developed an efficiency statistic based on bias and variance of the free energy estimates. For the two small OA binders, the free energy estimates computed with alchemical and potential of mean force approaches show relatively similar variance and bias as a function of the number of energy/force evaluations, with the attach-pull-release (APR), GROMACS expanded ensemble, and NAMD double decoupling submissions obtaining the greatest efficiency. The differences between the methods increase when analyzing the CB8-quinine system, where both the guest size and correlation times for system dynamics are greater. For this system, nonequilibrium switching (GROMACS/NS-DS/SB) obtained the overall highest efficiency. Surprisingly, the results suggest that specifying force field parameters and partial charges is insufficient to generally ensure reproducibility, and we observe differences between seemingly converged predictions ranging approximately from 0.3 to 1.0 kcal/mol, even with almost identical simulations parameters and system setup (e.g., Lennard-Jones cutoff, ionic composition). Further work will be required to completely identify the exact source of these discrepancies. Among the conclusions emerging from the data, we found that Hamiltonian replica exchange-while displaying very small variance-can be affected by a slowly-decaying bias that depends on the initial population of the replicas, that bidirectional estimators are significantly more efficient than unidirectional estimators for nonequilibrium free energy calculations for systems considered, and that the Berendsen barostat introduces non-negligible artifacts in expanded ensemble simulations.

Pyruvate carboxylase deficiency type A and type C: Characterization of five novel pathogenic variants in PC and analysis of the genotype-phenotype correlation.

Pyruvate carboxylase deficiency (PCD) is caused by biallelic mutations of the PC gene. The reported clinical spectrum includes a neonatal form with early death (type B), an infantile fatal form (type A), and a late-onset form with isolated mild intellectual delay (type C). Apart from homozygous stop-codon mutations leading to type B PCD, a genotype-phenotype correlation has not otherwise been discernible. Indeed, patients harboring biallelic heterozygous variants leading to PC activity near zero can present either with a fatal infantile type A or with a benign late onset type C form. In this study, we analyzed six novel patients with type A (three) and type C (three) PCD, and compared them with previously reported cases. First, we observed that type C PCD is not associated to homozygous variants in PC. In silico modeling was used to map former and novel variants associated to type A and C PCD, and to predict their potential effects on the enzyme structure and function. We found that variants lead to type A or type C phenotype based on the destabilization between the two major enzyme conformers. In general, our study on novel and previously reported patients improves the overall understanding on type A and C PCD.

Performance evaluation of molecular docking and free energy calculations protocols using the D3R Grand Challenge 4 dataset.

Using the D3R Grand Challenge 4 dataset containing Beta-secretase 1 (BACE) and Cathepsin S (CatS) inhibitors, we have evaluated the performance of our in-house docking workflow that involves in the first step the selection of the most suitable docking software for the system of interest based on structural and functional information available in public databases, followed by the docking of the dataset to predict the binding modes and ranking of ligands. The macrocyclic nature of the BACE ligands brought additional challenges, which were dealt with by a careful preparation of the three-dimensional input structures for ligands. This provided top-performing predictions for BACE, in contrast with CatS, where the predictions in the absence of guiding constraints provided poor results. These results highlight the importance of previous structural knowledge that is needed for correct predictions on some challenging targets. After the end of the challenge, we also carried out free energy calculations (i.e. in a non-blinded manner) for CatS using the pmx software and several force fields (AMBER, Charmm). Using knowledge-based starting pose construction allowed reaching remarkable accuracy for the CatS free energy estimates. Interestingly, we show that the use of a consensus result, by averaging the results from different force fields, increases the prediction accuracy.

Utilizing dipole-dipole cross-correlated relaxation for the measurement of angles between pairs of opposing CαHα-CαHα bonds in anti-parallel ß-sheets.

Dipole-dipole cross-correlated relaxation (CCR) between two spin pairs is rich with macromolecular structural and dynamic information on inter-nuclear bond vectors. Measurement of short range dipolar CCR rates has been demonstrated for a variety of inter-nuclear vector spin pairs in proteins and nucleic acids, where the multiple quantum coherence necessary for observing the CCR rate is created by through-bond scalar coupling. In principle, CCR rates can be measured for any pair of inter-nuclear vectors where coherence can be generated between one spin of each spin pair, regardless of both the distance between the two spin pairs and the distance of the two spins forming the multiple quantum coherence. In practice, however, long range CCR (lrCCR) rates are challenging to measure due to difficulties in linking spatially distant spin pairs. By utilizing through-space relaxation allowed coherence transfer (RACT), we have developed a new method for the measurement of lrCCR rates involving CαHα bonds on opposing anti-parallel ß-strands. The resulting lrCCR rates are straightforward to interpret since only the angle between the two vectors modulates the strength of the interference effect. We applied our lrCCR measurement to the third immunoglobulin-binding domain of the streptococcal protein G (GB3) and utilize published NMR ensembles and static NMR/X-ray structures to highlight the relationship between the lrCCR rates and the CαHα-CαHα inter-bond angle and bond mobility. Furthermore, we employ the lrCCR rates to guide the selection of sub-ensembles from the published NMR ensembles for enhancing the structural and dynamic interpretation of the data. We foresee this methodology for measuring lrCCR rates as improving the generation of structural ensembles by providing highly accurate details concerning the orientation of CαHα bonds on opposing anti-parallel ß-strands.

Insights into the molecular basis for substrate binding and specificity of the wild-type L-arginine/agmatine antiporter AdiC.

Pathogenic enterobacteria need to survive the extreme acidity of the stomach to successfully colonize the human gut. Enteric bacteria circumvent the gastric acid barrier by activating extreme acid-resistance responses, such as the arginine-dependent acid resistance system. In this response, l-arginine is decarboxylated to agmatine, thereby consuming one proton from the cytoplasm. In Escherichia coli, the l-arginine/agmatine antiporter AdiC facilitates the export of agmatine in exchange of l-arginine, thus providing substrates for further removal of protons from the cytoplasm and balancing the intracellular pH. We have solved the crystal structures of wild-type AdiC in the presence and absence of the substrate agmatine at 2.6-Å and 2.2-Å resolution, respectively. The high-resolution structures made possible the identification of crucial water molecules in the substrate-binding sites, unveiling their functional roles for agmatine release and structure stabilization, which was further corroborated by molecular dynamics simulations. Structural analysis combined with site-directed mutagenesis and the scintillation proximity radioligand binding assay improved our understanding of substrate binding and specificity of the wild-type l-arginine/agmatine antiporter AdiC. Finally, we present a potential mechanism for conformational changes of the AdiC transport cycle involved in the release of agmatine into the periplasmic space of E. coli.

pmx Webserver: A User Friendly Interface for Alchemistry.

With the increase of available computational power and improvements in simulation algorithms, alchemical molecular dynamics based free energy calculations have developed into routine usage. To further facilitate the usability of alchemical methods for amino acid mutations, we have developed a web based infrastructure for obtaining hybrid protein structures and topologies. The presented webserver allows amino acid mutation selection in five contemporary molecular mechanics force fields. In addition, a complete mutation scan with a user defined amino acid is supported. The output generated by the webserver is directly compatible with the Gromacs molecular dynamics engine and can be used with any of the alchemical free energy calculation setup. Furthermore, we present a database of input files and precalculated free energy differences for tripeptides approximating a disordered state of a protein, of particular use for protein stability studies. Finally, the usage of the webserver and its output is exemplified by performing an alanine scan and investigating thermodynamic stability of the Trp cage mini protein. The webserver is accessible at http://pmx.mpibpc.mpg.de.

Accurate and Rigorous Prediction of the Changes in Protein Free Energies in a Large-Scale Mutation Scan.

The prediction of mutation-induced free-energy changes in protein thermostability or protein-protein binding is of particular interest in the fields of protein design, biotechnology, and bioengineering. Herein, we achieve remarkable accuracy in a scan of 762 mutations estimating changes in protein thermostability based on the first principles of statistical mechanics. The remaining error in the free-energy estimates appears to be due to three sources in approximately equal parts, namely sampling, force-field inaccuracies, and experimental uncertainty. We propose a consensus force-field approach, which, together with an increased sampling time, leads to a free-energy prediction accuracy that matches those reached in experiments. This versatile approach enables accurate free-energy estimates for diverse proteins, including the prediction of changes in the melting temperature of the membrane protein neurotensin receptor 1.

Improved validation of IDP ensembles by one-bond Cα-Hα scalar couplings.

Intrinsically disordered proteins (IDPs) are best described by ensembles of conformations and a variety of approaches have been developed to determine IDP ensembles. Because of the large number of conformations, however, cross-validation of the determined ensembles by independent experimental data is crucial. The (1)JCαHα coupling constant is particularly suited for cross-validation, because it has a large magnitude and mostly depends on the often less accessible dihedral angle ψ. Here, we reinvestigated the connection between (1)JCαHα values and protein backbone dihedral angles. We show that accurate amino-acid specific random coil values of the (1)JCαHα coupling constant, in combination with a reparameterized empirical Karplus-type equation, allow for reliable cross-validation of molecular ensembles of IDPs.

pmx: Automated protein structure and topology generation for alchemical perturbations.

Computational protein design requires methods to accurately estimate free energy changes in protein stability or binding upon an amino acid mutation. From the different approaches available, molecular dynamics-based alchemical free energy calculations are unique in their accuracy and solid theoretical basis. The challenge in using these methods lies in the need to generate hybrid structures and topologies representing two physical states of a system. A custom made hybrid topology may prove useful for a particular mutation of interest, however, a high throughput mutation analysis calls for a more general approach. In this work, we present an automated procedure to generate hybrid structures and topologies for the amino acid mutations in all commonly used force fields. The described software is compatible with the Gromacs simulation package. The mutation libraries are readily supported for five force fields, namely Amber99SB, Amber99SB*-ILDN, OPLS-AA/L, Charmm22*, and Charmm36.

Guidelines for Free-Energy Calculations Involving Charge Changes.

The Coulomb interactions in molecular simulations are inherently approximated due to the finite size of the molecular box sizes amenable to current-day compute power. Several methods exist for treating long-range electrostatic interactions, yet these approaches are subject to various finite-size-related artifacts. Lattice-sum methods are frequently used to approximate long-range interactions; however, these approaches also suffer from artifacts which become particularly pronounced for free-energy calculations that involve charge changes. The artifacts, however, also affect the sampling when plain simulations are performed, leading to a biased ensemble. Here, we investigate two previously described model systems to determine if artifacts continue to play a role when overall neutral boxes are considered, in the context of both free-energy calculations and sampling. We find that ensuring that no net-charge changes take place, while maintaining a neutral simulation box, may be sufficient provided that the simulation boxes are large enough. Addition of salt to the solution (when appropriate) can further alleviate the remaining artifacts in the sampling or the calculated free-energy differences. We provide practical guidelines to avoid finite-size artifacts.

Optimal superpositioning of flexible molecule ensembles.

Analysis of the internal dynamics of a biological molecule requires the successful removal of overall translation and rotation. Particularly for flexible or intrinsically disordered peptides, this is a challenging task due to the absence of a well-defined reference structure that could be used for superpositioning. In this work, we started the analysis with a widely known formulation of an objective for the problem of superimposing a set of multiple molecules as variance minimization over an ensemble. A negative effect of this superpositioning method is the introduction of ambiguous rotations, where different rotation matrices may be applied to structurally similar molecules. We developed two algorithms to resolve the suboptimal rotations. The first approach minimizes the variance together with the distance of a structure to a preceding molecule in the ensemble. The second algorithm seeks for minimal variance together with the distance to the nearest neighbors of each structure. The newly developed methods were applied to molecular-dynamics trajectories and normal-mode ensembles of the Aß peptide, RS peptide, and lysozyme. These new (to our knowledge) superpositioning methods combine the benefits of variance and distance between nearest-neighbor(s) minimization, providing a solution for the analysis of intrinsic motions of flexible molecules and resolving ambiguous rotations.

Computational analysis of local membrane properties.

In the field of biomolecular simulations, dynamics of phospholipid membranes is of special interest. A number of proteins, including channels, transporters, receptors and short peptides are embedded in lipid bilayers and tightly interact with phospholipids. While the experimental measurements report on the spatial and/or temporal average membrane properties, simulation results are not restricted to the average properties. In the current study, we present a collection of methods for an efficient local membrane property calculation, comprising bilayer thickness, area per lipid, deuterium order parameters, Gaussian and mean curvature. The local membrane property calculation allows for a direct mapping of the membrane features, which subsequently can be used for further analysis and visualization of the processes of interest. The main features of the described methods are highlighted in a number of membrane systems, namely: a pure dimyristoyl-phosphatidyl-choline (DMPC) bilayer, a fusion peptide interacting with a membrane, voltage-dependent anion channel protein embedded in a DMPC bilayer, cholesterol enriched bilayer and a coarse grained simulation of a curved palmitoyl-oleoyl-phosphatidyl-choline lipid membrane. The local membrane property analysis proves to provide an intuitive and detailed view on the observables that are otherwise interpreted as averaged bilayer properties.

Accurately Predicting Protein pKa Values Using Nonequilibrium Alchemy.

The stability, solubility, and function of a protein depend on both its net charge and the protonation states of its individual residues. pKa is a measure of the tendency for a given residue to (de)protonate at a specific pH. Although pKa values can be resolved experimentally, theory and computation provide a compelling alternative. To this end, we assess the applicability of a nonequilibrium (NEQ) alchemical free energy method to the problem of pKa prediction. On a data set of 144 residues that span 13 proteins, we report an average unsigned error of 0.77 ± 0.09, 0.69 ± 0.09, and 0.52 ± 0.04 pK for aspartate, glutamate, and lysine, respectively. This is comparable to current state-of-the-art predictors and the accuracy recently reached using free energy perturbation methods (e.g., FEP+). Moreover, we demonstrate that our open-source, pmx-based approach can accurately resolve the pKa values of coupled residues and observe a substantial performance disparity associated with the lysine partial charges in Amber14SB/Amber99SB*-ILDN, for which an underused fix already exists.