NE-RDFE: A protocol and toolkit for computing relative dissociation free energies with GROMACS between dissimilar molecules using bidirectional nonequilibrium dual topology schemes.

We describe a step-by-step protocol and toolkit for the computation of the relative dissociation free energy (RDFE) with the GROMACS molecular dynamics package, based on a novel bidirectional nonequilibrium alchemical approach. The proposed methodology does not require any intervention on the code and allows computing with good accuracy the RDFE between small molecules with arbitrary differences in volume, charge, and chemical topology. The procedure is illustrated for the challenging SAMPL9 batch of host-guest pairs. The article is supplemented by a detailed online tutorial, available at https://procacci.github.io/vdssb_gromacs/NE-RDFE and by a public Zenodo repository available at https://zenodo.org/record/6982932.

SAMPL9 blind predictions for toluene/water partition coefficients using nonequilibrium alchemical approaches.

We present our blind prediction of the toluene-water partition coefficients in the context of the SAMPL9 challenge. For the calculation of the solvation free energies in water, toluene, and 1-octanol, we used an efficient MD-based nonequilibrium alchemical technique relying on the GAFF2 non-polarizable force field. The method is based on the fast-growth of an initially decoupled solute. Canonical sampling of the associated end-state is efficiently obtained by performing a Hamiltonian replica exchange simulation of the gas-phase solute molecule alone, combined with equilibrium configurations of the solvent. Before submitting the prediction, a pre-assessment of the method and of the force field was made by comparing with the known experimental counterpart the calculated octanol-water partition coefficients using different set of atomic charges. The analysis allowed to optimize our blind prediction for the toluene-water partition coefficients, providing at the same time valid clues for improving the performance and reliability of the non-polarizable force field in free energy calculations of drug-receptor systems.

2-Butanol Aqueous Solutions: A Combined Molecular Dynamics and Small/Wide-Angle X-ray Scattering Study.

Structural properties of 2-butanol aqueous solutions at different concentrations have been studied using small- and wide-angle X-ray scattering and molecular dynamics simulations. The experimental structure factors have been accurately reproduced by the simulations, allowing one to explain their variation with concentration and to achieve a detailed description of the structural and dynamic properties of the studied systems. The analysis of experimental and computational data has shown that 2-butanol, the simplest aliphatic chiral alcohol, tends to form aggregates at a concentration above 1 M, affecting also both the structural and dynamic properties of the solvent.

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.

Does Hamiltonian Replica Exchange via Lambda-Hopping Enhance the Sampling in Alchemical Free Energy Calculations?

In the context of computational drug design, we examine the effectiveness of the enhanced sampling techniques in state-of-the-art free energy calculations based on alchemical molecular dynamics simulations. In a paradigmatic molecule with competition between conformationally restrained E and Z isomers whose probability ratio is strongly affected by the coupling with the environment, we compare the so-called λ-hopping technique to the Hamiltonian replica exchange methods assessing their convergence behavior as a function of the enhanced sampling protocols (number of replicas, scaling factors, simulation times). We found that the pure λ-hopping, commonly used in solvation and binding free energy calculations via alchemical free energy perturbation techniques, is ineffective in enhancing the sampling of the isomeric states, exhibiting a pathological dependence on the initial conditions. Correct sampling can be restored in λ-hopping simulation by the addition of a "hot-zone" scaling factor to the λ-stratification (FEP+ approach), provided that the additive hot-zone scaling factors are tuned and optimized using preliminary ordinary replica-exchange simulation of the end-states.

Virtual Double-System Single-Box for Absolute Dissociation Free Energy Calculations in GROMACS.

We describe a step-by-step protocol for the computation of absolute dissociation free energy with GROMACS code and PLUMED library, which exploits a combination of advanced sampling techniques and nonequilibrium alchemical methodologies. The computational protocol has been automated through an open source Python middleware (HPC_Drug) which allows one to set up the GROMACS/PLUMED input files for execution on high performing computing facilities. The proposed protocol, by exploiting its inherent parallelism and the power of the GROMACS code on graphical processing units, has the potential to afford accurate and precise estimates of the dissociation constants in drug-receptor systems described at the atomistic level. The procedure has been applied to the calculation of the absolute dissociation free energy of PF-07321332, an oral antiviral proposed by Pfizer, with the main protease (3CLpro) of SARS-CoV-2.

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 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.

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.

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 blind predictions of water-octanol partition coefficients using nonequilibrium alchemical approaches.

In this paper, we compute, by means of a non equilibrium alchemical technique, the water-octanol partition coefficients (LogP) for a series of drug-like compounds in the context of the SAMPL6 challenge initiative. Our blind predictions are based on three of the most popular non-polarizable force fields, CGenFF, GAFF2, and OPLS-AA and are critically compared to other MD-based predictions produced using free energy perturbation or thermodynamic integration approaches with stratification. The proposed non-equilibrium method emerges has a reliable tool for LogP prediction, systematically being among the top performing submissions in all force field classes for at least two among the various indicators such as the Pearson or the Kendall correlation coefficients or the mean unsigned error. Contrarily to the widespread equilibrium approaches, that yielded apparently very disparate results in the SAMPL6 challenge, all our independent prediction sets, irrespective of the adopted force field and of the adopted estimate (unidirectional or bidirectional) are, mutually, from moderately to strongly correlated.

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 3CLpro 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.

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 sp3 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.

Precision and computational efficiency of nonequilibrium alchemical methods for computing free energies of solvation. II. Unidirectional estimates.

The present paper is the second part of a series of papers aimed at assessing the accuracy of alchemical computational approaches based on nonequilibrium techniques for solvation free energy of organic molecules in the context of molecular dynamics simulations. In Paper I [Procacci, J. Chem. Phys. 151, 144113 (2019)], we dealt with bidirectional estimates of solvation free energies using nonequilibrium approaches. Here, we assess accuracy and precision of unidirectional estimates with the focus on the Gaussian and Jarzynski estimators. We present a very simple methodology to increase the statistics in the work distribution, hence boosting the accuracy and precision of the Jarzynski unidirectional estimates at no extra cost, exploiting the observed decorrelation between the random variables represented by the Lennard-Jones solute-solvent recoupling or decoupling work and by the electrostatic work due to the charging/discharging of the solute in the solvent.

Accuracy, precision, and efficiency of nonequilibrium alchemical methods for computing free energies of solvation. I. Bidirectional approaches.

In the context of molecular dynamics simulations, alchemical approaches based on nonequilibrium techniques are recently emerging as a powerful method for the computation of solvation free energy of druglike compounds. Here, we present a rigorous and extensive analysis of the accuracy and precision of the method as a function of the parameters qualifying the nonequilibrium alchemical protocol (e.g., number and length of the nonequilibrium trajectories and decoupling or recoupling alchemical schedule) on a selection of drug-size organic compounds characterized by a nontrivial conformational free energy landscape. The study is organized in two contributions. The first paper includes a detailed description of method and of the conformational behavior of molecular systems. Results are focused on the accuracy and precision bidirectional estimates of solvation free energy, notably those based on the so-called Bennett acceptance ratio. In the second paper, unidirectional estimates for solvation free energy are analyzed in depth.

Blocking the FKBP12 induced dendrimeric burst in aberrant aggregation of α-synuclein by using the ElteN378 synthetic inhibitor.

α-Synuclein (α-syn), a disordered cytoplasmatic protein, plays a fundamental role in the pathogenesis of Parkinson's disease (PD). Here, we have shown, using photophysical measurements, that addition of FKBP12 to α-syn solutions, dramatically accelerates protein aggregation, leading to an explosion of dendritic structures revealed by fluorescence and phase-contrast microscopy. We have further demonstrated that this aberrant α-syn aggregation can be blocked using a recently discovered non-immunosuppressive synthetic inhibitor of FKBP12, ElteN378. The role of FKBP12 and of ElteN378 in the α-syn aggregation mechanism has been elucidated using molecular dynamics simulations based on an effective coarse-grained model. The reported data not only reveal a new potent synthetic drug as a candidate for early stage treatment of α-syn dependent neurodegenerations but also pave the way to a deeper understanding of the mechanism of action of FKBP12 on α-syn oligomeric aggregation, a topic which is still controversial.

Lipid tempering simulation of model biological membranes on parallel platforms.

In this report we have tested a parallel implementation for the simulation of lipid bilayers at the atomistic level, based on a generalized ensemble protocol where only the torsional degrees of freedom of the alkyl chains of the lipids are heated. The results in terms of configurational sampling enhancement have been compared with a conventional simulation produced with a widespread molecular dynamics code. Results show that the proposed thermodynamic-based multiple trajectories parallel protocol for membrane simulations allows for an efficient use of CPU resources with respect to the conventional single trajectory, providing accurate results for area and volume per lipid, membrane thickness, undulation spectra and boosting significantly diffusion and mixing in lipid bilayers due to the sampling enhancement of gauche/trans ratios of the alkyl chain dihedral angles.

SAMPL6 host-guest blind predictions using a non equilibrium alchemical approach.

In this paper, we compute, by means of a non equilibrium alchemical technique, called fast switching double annihilation methods (FSDAM), the absolute standard dissociation free energies of the the octa acids host-guest systems in the SAMPL6 challenge initiative. FSDAM is based on the production of canonical configurations of the bound and unbound states via enhanced sampling and on the subsequent generation of hundreds of fast non-equilibrium ligand annihilation trajectories. The annihilation free energies of the ligand when bound to the receptor and in bulk solvent are obtained from the collection of work values using an estimate based on the Crooks theorem for driven non equilibrium processes. The FSDAM blind prediction, relying on the normality assumption for the annihilation work distributions, ranked fairly well among the submitted blind predictions that were not adjusted with a linear corrections obtained from retrospective data on similar host guest systems. Improved results for FSDAM can be obtained by post-processing the work data assuming mixtures of normal components.

Correspondence between light-absorption spectrum and nonequilibrium work distribution as a mean to access free energy differences between electronic states.

The problem of recovering the free energy difference between two electronic states has been investigated by Frezzato [Chem. Phys. Lett. 533, 106 (2012)], exploring the equivalence between light-absorption spectra and work distribution, hence opening to the application of a spectroscopic version of the Jarzynski equality (JE) [Phys. Rev. Lett. 78, 2690 (1997)]. Here, assuming the validity of the time-dependent perturbation theory, we demonstrate that such equivalence does not lead to the known form of the JE. This is ascribed to the fact that light-absorption processes cannot be described as stochastic processes. To emphasize such an aspect, we devise a stochastic model for the UV-vis (ultraviolet and visible) absorption, suitable for determining the free energy difference between two generic quantum manifolds in a JE-like fashion. However, the model would require explicit knowledge of the transition dipole moments, which are in general not available. Nonetheless, we derive a spectroscopic version of the JE that allows us to recover the free energy difference between the ground and an excited electronic state when the latter state is the only one observed in the spectrum.