Search

In the companion article (Giovannelli et al., 10.1021/acs.jctc.7b00594), we presented an alchemical approach, based on nonequilibrium molecular dynamics simulations, to compute absolute binding free energies of a generic host-guest system. Two alternative computational routes, called binded-domain and single-point alchemical-path schemes, have been proposed. This study is addressed to furnish numerical validation and illustrative examples of the above-mentioned alchemical schemes. Validation is provided by comparing binding free-energy data relative to two poses of a Zn(II)·anion complex with those recovered from an alternative approach, based on steered molecular dynamics simulations. We illustrate important technical and theoretical aspects for a good practice in applying both alchemical schemes, not only through the calculations on the Zn(II)·anion complex, but also estimating absolute binding free energies of 1:1 complexes of β-cyclodextrin with aromatic compounds (benzene and naphthalene). Comparison with experimental data and previous molecular dynamics simulation studies further confirms the validity of the present nonequilibrium-alchemical methodology.

The fast-switching decoupling method is a powerful nonequilibrium technique to compute absolute binding free energies of ligand-receptor complexes (Sandberg et al., J. Chem. Theory Comput. 2014, 11, 423-435). Inspired by the theory of noncovalent binding association of Gilson and co-workers (Biophys. J. 1997, 72, 1047-1069), we develop two approaches, termed binded-domain and single-point alchemical-path schemes (BiD-AP and SiP-AP), based on the possibility of performing alchemical trajectories during which the ligand is constrained to fixed positions relative to the receptor. The BiD-AP scheme exploits a recent generalization of nonequilibrium work theorems to estimate the free energy difference between the coupled and uncoupled states of the ligand-receptor complex. With respect to the fast-switching decoupling method without constraints, BiD-AP prevents the ligand from leaving the binding site, but still requires an estimate of the positional binding-site volume, which may not be a simple task. On the other side, the SiP-AP scheme allows avoidance of the calculation of the binding-site volume by introducing an additional equilibrium simulation of ligand and receptor in the bound state. In the companion article (DOI: 10.1021/acs.jctc.7b00595), we show that the extra computational effort required by SiP-AP leads to a significant improvement of accuracy in the free energy estimates.

An important issue concerning computer simulations addressed to free energy estimates via nonequilibrium work theorems, such as the Jarzynski equality [Phys. Rev. Lett. 1997, 78, 2690], is the computational effort required to achieve results with acceptable accuracy. In this respect, the dynamical freezing approach [Phys. Rev. E 2009, 80, 041124] has been shown to improve the efficiency of this kind of simulations, by blocking the dynamics of particles located outside an established mobility region. In this report, we show that dynamical freezing produces a systematic spurious decrease of the particle density inside the mobility region. As a consequence, the requirements to apply nonequilibrium work theorems are only approximately met. Starting from these considerations, we have developed a simulation scheme, called "elastic barrier dynamical freezing", according to which a stiff potential-energy barrier is enforced at the boundaries of the mobility region, preventing the particles from leaving this region of space during the nonequilibrium trajectories. The method, tested on the calculation of the distance-dependent free energy of a dimer immersed into a Lennard-Jones fluid, provides an accuracy comparable to the conventional steered molecular dynamics, with a computational speedup exceeding a few orders of magnitude.

A simulation-based approach is proposed to estimate free energy differences between configurational states A and B, defined in terms of collective coordinates of the molecular system. The computational protocol is organized into three stages that can be carried on simultaneously. Two of them consist of independent simulations aimed at sampling, in turn, A and B states. In order to limit the evolution of the system around A and B, biased sampling simulations such as umbrella sampling can be employed. These simulations allow us to estimate local configuration integrals associated with A and B, which can be viewed as vibrational contributions to the free energy. Free energy evaluation is completed by the linking-path stage, in which the potential of mean force difference is estimated between two arbitrary points of the configurational surface, located the first around A and the second around B. The linking path in the space of the collective coordinates is arbitrary and can be computed with any method, starting from adaptive biasing potential/force approaches to nonequilibrium techniques. As an illustrative example, we present the calculation of free energy differences between conformational states of the alanine dipeptide in the space of backbone dihedral angles. The basic advantage of this method, that we term "path-linked domains" scheme, is to prevent accurate calculation of the whole free energy hypersurface in the space of the collective coordinates, thus limiting the statistical sampling to a minimum. Path-linked domains schemes can be applied to a variety of biochemical processes, such as protein-ligand complexation or folding-unfolding interconversion.

Nonequilibrium Candidate Monte Carlo simulation [Nilmeier et al., Proc. Natl. Acad. Sci. U.S.A. 2011, 108, E1009-E1018] is a tool devised to design Monte Carlo moves with high acceptance probabilities that connect uncorrelated configurations. Such moves are generated through nonequilibrium driven dynamics, producing candidate configurations accepted with a Monte Carlo-like criterion that preserves the equilibrium distribution. The probability of accepting a candidate configuration as the next sample in the Markov chain basically depends on the work performed on the system during the nonequilibrium trajectory and increases with decreasing such a work. It is thus strategically relevant to find ways of producing nonequilibrium moves with low work, namely moves where dissipation is as low as possible. This is the goal of our methodology, in which we combine Nonequilibrium Candidate Monte Carlo with Configurational Freezing schemes developed by Nicolini et al. (J. Chem. Theory Comput. 2011, 7, 582-593). The idea is to limit the configurational sampling to particles of a well-established region of the simulation sample, namely the region where dissipation occurs, while leaving fixed the other particles. This allows to make the system relaxation faster around the region perturbed by the finite-time switching move and hence to reduce the dissipated work, eventually enhancing the probability of accepting the generated move. Our combined approach enhances significantly configurational sampling, as shown by the case of a bistable dimer immersed in a dense fluid.

Wavelet transform has been used to correlate spectroscopic and structural properties from trajectories obtained by ab initio molecular dynamics simulations. This method has been applied to hydrogen bond dynamics of glycols in heavy water solutions, showing how the stretching frequency of the intramolecular O-H bond changes with the intermolecular hydrogen-bond distance. The resulting wavelet spectrograms have been interpreted according to H-bond strength and stability.

The reaction of ethylene in condensed phases under high pressure has been investigated by ab initio molecular dynamics. Both disordered and crystalline samples have been simulated, and some insights on the reaction mechanism have been obtained. System size effects have been investigated for the disordered samples. A polymerization reaction occurs by an ionic mechanism. In both the disordered and the crystal phases, the reaction products obtained (linear chains in the disordered systems and branched chains in the crystal) are in qualitative agreement with the experiments.

Ab initio Car-Parrinello molecular dynamics simulations have been performed in order to investigate the solvation properties of Mg(2+) and Ca(2+) in fully deuterated methanol solution to better understand polarization effects induced by the ions. Charge transfer and dipole moment calculations have been performed to give more detailed insight on the role of the electronic reorganization and its effect on the first solvation shell stability. The perturbation of the methanol H-bond network has been investigated.