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.

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

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.

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.

Structural properties of 2-butanol aqueous solutions at different concentrations have been studied using small-and wide-angle Xray 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.

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.

This work proposes a voltammetric aptasensor to detect deoxynivalenol (DON) mycotoxin. The development steps of the aptasensor were partnered for the first time to a computational study to gain insights onto the molecular mechanisms involved into the interaction between a thiol-tethered DNA aptamer (80mer-SH) and DON. The exploited docking study allowed to find the binding region of the oligonucleotide sequence and to determine DON preferred orientation. A biotinylated oligonucleotide sequence (20mer-BIO) complementary to the aptamer was chosen to carry out a competitive format. Graphite screen-printed electrodes (GSPEs) were electrochemically modified with polyaniline and gold nanoparticles (AuNPs@PANI) by means of cyclic voltammetry (CV) and worked as a scaffold for the immobilization of the DNA aptamer. Solutions containing increasing concentrations of DON and a fixed amount of 20mer-BIO were dropped onto the aptasensor surface: the resulting hybrids were labeled with an alkaline phosphatase (ALP) conjugate to hydrolyze 1-naphthyl phosphate (1-NPP) substrate into 1-naphthol product, detected by differential pulse voltammetry (DPV). According to its competitive format, the aptasensor response was signal-off in the range 5.0--30.0 ngmL?1 DON. A detection limit of 3.2 ngmL?1 was achieved within a 1-hour detection time. Preliminary experiments on maize flour samples spiked with DON yielded good recovery values.

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.

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.

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.

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.

Using a combination of enhanced sampling molecular dynamics techniques and non-equilibrium alchemical transformations with full atomistic details, we have shown that hydroxychloroquine (HCQ) may act as a mild inhibitor of important functional proteins for SARS-CoV2 replication, with potency increasing in the series PLpro, 3CLpro, RdRp. By analyzing the bound state configurations, we were able to improve the potency for the 3CLpro target, designing a novel HCQ-inspired compound, named PMP329, with predicted nanomolar activity. If confirmed in vitro, our results provide a molecular rationale for the use of HCQ or of strictly related derivatives in the treatment of Covid-19. This journal is

ABSTRACT: 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 (3CLpro) of SARS-CoV-2. The core structures of 3CLpro 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 3CLpro 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 3CLpro inhibition, which may constitute possible lead compounds for the development of antiviral agents in Covid-19 treatment.

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 kBT, with the potential of producing artifacts and entanglements if not judiciously implemented.

Intrinsically disordered proteins, such as tau or α-synuclein, have long been associated with a dysfunctional role in neurodegenerative diseases. In Alzheimer's and Parkinson's' diseases, these proteins, sharing a common chemical-physical pattern with alternating hydrophobic and hydrophilic domains rich in prolines, abnormally aggregate in tangles in the brain leading to progressive loss of neurons. In this review, we present an overview linking the studies on the implication of the peptidyl-prolyl isomerase domain of immunophilins, and notably FKBP12, to a variety of neurodegenerative diseases, focusing on the molecular origin of such a role. The involvement of FKBP12 dysregulation in the aberrant aggregation of disordered proteins pinpoints this protein as a possible therapeutic target and, at the same time, as a predictive biomarker for early diagnosis in neurodegeneration, calling for the development of reliable, fast and cost-effective detection methods in body fluids for community-based screening campaigns.

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.

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.;et al. 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.

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.

Water has a fundamental role in important processes spanning a wide range of pressure and temperature conditions. Knowledge of structural, dynamic and thermodynamic properties of water at nonstandard conditions is a primary concern since interest in astronomical, geological, and technological processes is continuously growing. Molecular dynamics simulations allow us to study thermodynamic conditions that require sophisticated techniques and instruments, while at the same time offering the interpretation of properties at the atomic level. It is established that the behavior of water is strongly affected by the temperature and pressure conditions, determining the existence of low and high density regimes. For the first time, a thermodynamic property, isothermal compressibility, has been adopted to detect the low-high density turning point at ambient temperature in liquid water due to pressure. Molecular dynamics simulations have been performed with five three-site models, allowing us to characterize the complexity of water nature at these conditions at the atomic level.

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.

In this paper, we compute, by means of a non equilibrium alchemical technique, the water-octanol partition coefcients (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 felds, CGenFF, GAFF2, and OPLS-AA and are critically compared to other MD-based predictions produced using free energy perturbation or thermodynamic integration approaches with stratifcation. The proposed non-equilibrium method emerges has a reliable tool for LogP prediction, systematically being among the top performing submissions in all force feld classes for at least two among the various indicators such as the Pearson or the Kendall correlation coefcients 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 feld and of the adopted estimate (unidirectional or bidirectional) are, mutually, from moderately to strongly correlated.

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.

Molecular dynamics simulations have been performed to compute the solvation free energy and the octanol/water partition coefficients for a challenging set of selected organic molecules, characterized by the simultaneous presence of functional groups coarsely spanning a large portion of the chemical space in drug-like compounds and, in many cases, by a complex conformational landscape (2-propoxyethanol, acetylsalicylic acid, cyclohexanamine, dialifor, ketoprofen, nitralin, profluralin, terbacil). OPLS-AA and GAFF2 parameterizations of the organic molecules and of 1-octanol have been done via the web-based automatic parameter generators, LigParGen [Dodda et al. Nucl. Acids Res. 2017; 121, 3864] and PrimaDORAC [Procacci, J. Chem. Inf. Model. 2017; 57, 1240], respectively. For the water solvent, three popular three-point sites models, TIP3P, SPCE and OPC3, were tested. Solvation free energies in water and 1-octanol are evaluated using a recently developed non equilibrium alchemical technology [Procacci et al. J. Chem. Theory Comput. 2014; 10, 2813]. Extensive and accurate simulations including all possible combinations of organic molecule, solvent and solvent model, allowed to assess the accuracy with regard to solvation free energies of the latest release of two widespread force fields, OPLS and GAFF. The collected data are relevant in the evaluation of the predictive power of these classical force fields (and of the related support software for automated parameterization) with regard to binding free energies in drug-receptor system for industrial applications.

Comment on “Statistical efficiency of methods for computing free energy of hydration”[ J. Chem. Phys. 149, 144111 (2018)

The structural and dynamic properties of imidazole in aqueous solution have been studied by means of classical and ab initio molecular dynamics simulations. We developed a new force field for the imidazole molecule with improved modelling of the electrostatic interactions, specifically tailored to address the well known drawbacks of existing force fields based on the atomic fractional charges approach. To this end, we reparametrized the charge distribution on the heterocyclic ring, introducing an extra site accounting for the lone pair on the deprotonated nitrogen. The accuracy of the model in describing the hydrogen bond pattern in aqueous solvent has been confirmed by comparing the classical results on imidazole-water interactions to accurate Car-Parrinello molecular dynamics simulations. The proposed classical model reproduce satisfactorily the experimental water/octanol partition coefficient of imidazole, as well as the structure of the imidazole molecular crystal. The force field has been finally applied to simulate aqueous solutions at various imidazole concentrations to obtain information on both imidazole-water and imidazole-imidazole interactions, providing a description of the different molecular arrangements in solution.

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.

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

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.

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.

In this work, we compute, by means of a non-equilibrium alchemical technique (fast switching double annihilation methods, FSDAMs), the dissociation free energy for five recently discovered micromolar to sub-nanomolar inhibitors of the Myeloid cell leukemia 1 protein, a key regulator in cell survival and death, providing valuable clues in the chemical-physical determinants of Mcl-1 inhibition. Using the same methodology, we attempt the calculation of the dissociation free energy of the BH3 domain from PUMA protein, binding Mcl-1 in the α-helical state. The synthetic ligands have been parametrized using the recently released GAFF2 general force field [ http://ambermd.org ] by means of the automated assignment tool PrimaDORAC [ Procacci , P. J. Chem. Inf.

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.

Alchemical molecular dynamics simulations for evaluating the binding free energy in ligand-receptor systems are emerging a new powerful tool for in silico drug discovery projects. Nonetheless, theoretical and technical challenges for these methodologies remain high and their use in industrial applications is still limited. In this contribution, the many variants of the alchemical approach are critically revised, discussing their strengths as well as their pitfalls and entanglements and placing existing computational theories into the broader context of nonequilibrium thermodynamics.

We present PrimaDORAC, a simple and freely accessible web interface for generating the topology and the parameter files of organic or drug molecules, to be used in molecular mechanics or molecular dynamics calculations. The interface relies on our in-house FORTRAN90 parser, working on the recently released Generalized Amber Force Field parameter set (GAFF2). AM1/BCC charges are computed using the Public Domain MOPAC7 program and the bond charge corrections (BCC) reported in Ref. Jakalian et al, J Comp Chem, 2002, 23, 1623-1641. The interface has been tested on about 52000 compounds (identified with a CAS-registry number) taken from the National Cancer Institute (NCI) Open Database. PrimaDORAC has been found very reliable producing GAFF2 minimized structures bearing a mean Root Square Displacement of about 0.01-0.02 nm with respect to the original CORINA generated 3D NCI structures. As a demonstrative example, we release the full topology and parameter files, along with the AM1/BCC-GAFF2 computed <em> in vacuo </em> IR spectrum, for some recently discovered PARP/MCL1 inhibitors. The web interface and parser, including the sources, are part of the ORAC code [Procacci, J. Chem. Inf.

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.

The present paper is intended to be a comprehensive assessment and rationalization, from a statistical mechanics perspective, of existing alchemical theories for binding free energy calculations of ligand-receptor systems. In detail, the statistical mechanics foundation of noncovalent interactions in ligand-receptor systems is revisited, providing a unifying treatment that encompasses the most important variants in the alchemical approaches from the seminal double annihilation method [Jorgensen et al. J. Chem. Phys. 1988 ; 89 , 3742] to the double decoupling method [Gilson et al. Biophys. J. 1997 ; 72 , 1047] and the Deng and Roux alchemical theory [Deng and Roux J. Chem. Theory Comput. 2006; 2 , 1255]. Connections and differences between the various alchemical approaches are highlighted and discussed.

The recently proposed fast switching double annihilation (FS-DAM) [Cardelli et al., J. Chem. Theory Comput., 2015, 11, 423] is aimed at computing the absolute standard dissociation free energies for the chemical equilibrium RL ? R + L occurring in solution through molecular dynamics (MD) simulations at the atomistic level. The technique is based on the production of fast nonequilibrium annihilation trajectories of one of the species (the ligand) in the solvated RL complex and in the bulk solvent. As detailed in the companion theoretical paper, the free energies of these two nonequilibrium annihilation processes are recovered by using an unbiased unidirectional estimate derived from the Crooks theorem exploiting the inherent Gaussian nature of the annihilation work. The FS-DAM technique was successfully applied to the evaluation of the dissociation free energy of the complexes of Zn(II) cations with an inhibitor of the Tumor Necrosis Factor α converting enzyme. Here we apply the technique to a real drug--receptor system, by satisfactorily reproducing the experimental dissociation free energies of FK506-related bulky ligands towards the native FKBP12 enzyme and by predicting the dissociation constants for the same ligands towards the mutant I56D. The effect of such mutations on the binding affinity of FK506-related ligands is relevant for assessing the thermodynamic forces regulating molecular recognition in FKBP12 inhibition.

Correction to “Free Energy Reconstruction in Bidirectional Force Spectroscopy Experiments: The Effect of the Device Stiffness”

We have derived, in the context of the Rigid Rotor Harmonic Approximation (RRHO), a general mass and Planck's constant h independent expression for the dissociation free energy in ligand--receptor systems, featuring a systematically (anti-binding) additive negative entropic term depending on readily available ligand--receptor quantities. The proposed RRHO expression allows to straightforwardly compute the absolute standard dissociation free energy without resorting to expensive normal mode analysis or other dynamical matrix-based techniques for evaluating the entropic contribution, hence providing an effective scoring function for assessing docking poses with no adjustable parameters. Our RRHO formula was tested on a set of 55 ligand--receptor systems obtaining correlation coefficients and unsigned mean errors comparable to or better than those obtained with computationally demanding techniques for the dissociation entropy assessment. The proposed compact reformulation of the RRHO entropy term could constitute the basis for new and more effective scoring functions in molecular docking-based high-throughput virtual screening for drug discovery.

In this contribution I critically revise the alchemical reversible approach in the context of the statistical mechanics theory of non-covalent bonding in drug-receptor systems. I show that most of the pitfalls and entanglements for the binding free energy evaluation in computer simulations are rooted in the equilibrium assumption that is implicit in the reversible method. These critical issues can be resolved by using a non-equilibrium variant of the alchemical method in molecular dynamics simulations, relying on the production of many independent trajectories with a continuous dynamical evolution of an externally driven alchemical coordinate, completing the decoupling of the ligand in a matter of a few tens of picoseconds rather than nanoseconds. The absolute binding free energy can be recovered from the annihilation work distributions by applying an unbiased unidirectional free energy estimate, on the assumption that any observed work distribution is given by a mixture of normal distributions, whose components are identical in either direction of the non-equilibrium process, with weights regulated by the Crooks theorem. I finally show that the inherent reliability and accuracy of the unidirectional estimate of the decoupling free energies, based on the production of a few hundreds of non-equilibrium independent sub-nanosecond unrestrained alchemical annihilation processes, is a direct consequence of the funnel-like shape of the free energy surface in molecular recognition. An application of the technique to a real drug-receptor system is presented in the companion paper.

We present a new release (6.0β) of the ORAC program [Marsili et al. J. Comput. Chem. 2010, 31, 1106-1116] with a hybrid OpenMP/MPI (open multiprocessing message passing interface) multilevel parallelism tailored for generalized ensemble (GE) and fast switching double annihilation (FS-DAM) nonequilibrium technology aimed at evaluating the binding free energy in drug-receptor system on high performance computing platforms. The production of the GE or FS-DAM trajectories is handled using a weak scaling parallel approach on the MPI level only, while a strong scaling force decomposition scheme is implemented for intranode computations with shared memory access at the OpenMP level. The efficiency, simplicity, and inherent parallel nature of the ORAC implementation of the FS-DAM algorithm, project the code as a possible effective tool for a second generation high throughput virtual screening in drug discovery and design. The code, along with documentation, testing, and ancillary tools, is distributed under the provisions of the General Public License and can be freely downloaded at www.chim.unifi.it/orac .

The first example of one-electron oxidation of thia-bridged triarylamine heterohelicenes to the corresponding exceptionally stable radical cations, fully characterized, as hexafluoroantimonate salts, by means of UV-Vis, EPR, ENDOR, density functional theory calculations and X-ray analyses, is reported. Chemical and electrochemical reversible redox processes are solidly demonstrated

In this paper, we present an improved method for obtaining unbiased estimates of the free energy difference between two thermodynamic states using the work distribution measured in nonequilibrium driven experiments connecting these states. The method is based on the assumption that any observed work distribution is given by a mixture of Gaussian distributions, whose normal components are identical in either direction of the nonequilibrium process, with weights regulated by the Crooks theorem. Using the prototypical example for the driven unfolding/folding of deca-alanine, we show that the predicted behavior of the forward and reverse work distributions, assuming a combination of only two Gaussian components with Crooks derived weights, explains surprisingly well the striking asymmetry in the observed distributions at fast pulling speeds. The proposed methodology opens the way for a perfectly parallel implementation of Jarzynski-based free energy calculations in complex systems.

Small-molecule inhibitors of Tumor Necrosis Factor α Converting Enzyme (TACE) are a promising therapeutic tool for Rheumatoid Arthritis, Multiple Sclerosis and other autoimmune diseases. Here we report on an extensive chemical--physical analysis of the media effects in modulating the conformational landscape of MBET306, the common scaffold and a synthetic precursor of a family of recently discovered tartrate-based TACE inhibitors. The structural features of this molecule with potential pharmaceutical applications have been disclosed by interpreting extensive photophysical measurements in various solvents with the aid of enhanced sampling molecular dynamics simulations and time dependent density functional calculations. Using a combination of experimental and computational techniques, the paper provides a general protocol for studying the structure in solution of molecular systems characterized by the existence of conformational metastable states.

We introduce an effective technique for the calculation of the binding free energy in drug-receptor systems using nonequilibrium molecular dynamics and application of the Jarzynski theorem. In essence, this novel methodology constitutes the nonequilibrium adaptation of an ancient free energy perturbation technique called Double Annihilation Method, invented more than 25 years ago [J. Chem. Phys. 1988, 89, 3742--3746] and upon which modern approaches of binding free energy computation in drug-receptor systems are heavily based. The proposed computational approach, termed Fast Switching Double Annihilation Methods (FS-DAM) in honor of its ancient ancestor, is applied to a prototypical example system with multiple binding sites, proving its computational potential and versatility in unraveling multiple site or allosteric binding processes.

ORAC is a FORTRAN suite to simulate complex biosystems at the atomistic level. The program's engine is supplemented by multiple time steps integration, smooth particle mesh Ewald method, constant pressure and constant temperature algorithms. Quite recently, several advanced techniques for enhanced sampling in atomistic systems have been implemented, including replica exchange with solute tempering, metadynamics, steered molecular dynamics, serial generalized ensemble simulations and nonequilibrium alchemical transformations. All these computational technologies have been devised for parallel architectures using the standard MPI communication protocol.

The purpose of this study is to compare the gelling behavior of two molecules: a chiral compound and its achiral counterpart. The chiral partner is characterized by a rigid, chiral pyrrolidine nucleus, while the achiral one contains a ? exible diethanolamine moiety. The chiral compound is an already known good organogelator, but also the achiral compound shows remarkable gelling properties. Very interestingly, a small fraction of the chiral compound induces chirality and strong CD e ? ects in its aggregates with the achiral one. The observed chirality ampli ? cation corresponds to a peculiar sergeant-and-soldier e ? ect. Molecular modelling and CD calculations suggested a model for the supramolecular assembly of hetero- aggregates that ? ts the experimental data.

We present an efficient and rigorous implementation of the fast switching alchemical transformation for systems where electrostatic interactions are treated using the smooth particle mesh Ewald method. Free energies are computed using bidirectional nonequilibrium alchemical trajectories by applying the Crooks fluctuation theorem and the Bennett acceptance ratio to the collection of the final alchemical works. The technique is used for the evaluation of the 1-octanol/water partition coefficients for some selected organic molecules. Fast switching alchemical tranformations appear to be competitive, both in accuracy and in efficiency, with respect to the traditional methods based on independent equilibrium simulations of intermediate states.

We present a new computational strategy for calculating the absolute binding free energy for protein ligand association in the context of atomistic simulation in explicit solvent. The method is based on an appropriate definition of a solute tempering scheme enforced via Hamilton replica exchange method (HREM). The definition of “solute” includes both the ligand and the active site, with the remainder of the systems defined as “solvent”. The hydrophilicity of the solute and the solute torsional plus nonbonded intrasolute interactions are increased and decreased, respectively, along the replica progression, thus favoring the extrusion of the drug form the active site in the scaled states of the generalized ensemble. The proposed technique, named “Energy Driven Undocking” (EDU-HREM), completely bypasses the need for defining and/or identifying the relevant reaction coordinates in a ligand receptor interactions and allows the calculation of the absolute binding free energy in one single generalized simulation of the drug-receptor system. The methodology is applied, with encouraging results, to the calculation of the absolute binding free energy of some FK506-related ligands of the peptidyl prolyl cis-trans isomerase protein (FKBP12) with known dissociation constants. Aspects of the binding/inhibition mechanism in FKBP12 are also analyzed and discussed.

Due to its central role in immunosuppression and cell proliferation and due to its specific peptidyl-prolyl-isomerase (PPI) function, the FKBP protein family is at the crossroad of several important metabolic pathways. Members of this family, and notably FK506 binding protein (FKBP12), are thought to be involved in neurodegenerative diseases such as Alzheimer disease, Parkinson disease, multiple sclerosis, amyotrophic lateral sclerosis, as well as in proliferation disorders and cancer. Using an interdisciplinary approach based on computational, synthetic, and experimental techniques, we show that the best potential binders for FKBP proteins optimally expose the two contiguous carbonyl oxygen in the proline-mimetic chain for FKBP docking and are characterized by the abundance of rigid quasi-cyclic structures stabilized in aqueous solution by intraligand hydrophobic interactions mimicking the macrolide structure of the natural FKBP binders FK506 and Rapamycin. These peculiar structural and chemical-physical features define at the same time an ElteX compound and the minimal pharmacore in the FKBP family, shedding new light on the isomerization mechanism of the PPI domain. On the basis of the above hypothesis, we have successfully designed and synthesized several nanomolar ElteX FKBP12 ligands. Among these, ElteN378 is a new low atomic weight ligand with affinity comparable to that of the macrolide Rapamycin. 2013 American Chemical Society.

We have synthesized and done an extensive chemical--physical analysis of the behavior of a new compound, named MBET306, a synthetic precursor of the recently discovered tartrate-based inhibitors of the protein Tumor Necrosis factor-a Converting Enzyme (TACE). Experimental and theoretical data have shown that in water solution MBET306 is overwhelmingly found as a monoanion at physiological pH, in a conformation that differs substantially from that detected in the known co-crystal structures of MBET306 derivatives bound to TACE. The body of collected experimental and theoretical data indicates that the monoanionic species binds Zn(II) inducing a strong stabilization of the crystal-like arrangement of the central tartrate zinc-binding group, lending support for a two step TACE docking mechanism via a zinc-bound intermediate. The thorough chemical--physical characterization of the conformational behavior of free and zinc-bound MBET306 in water bulk solution opens new avenues for the rational drug design of tartrate-based highly specific TACE inhibitors

A two-step one purification access to dichloronaphtho[1,2-b:8,7-b] and [1,2-b:5,6-b]dithiophenes using bis-alkylnaphthyl alkynes and phthalimidesulfenyl chloride as starting materials has been developed. The functionalization of the carbon-chlorine bonds allowed further modification of NDT core, broadening the potential of the methodology. 2013 American Chemical Society.

In this paper, we illustrate a practical technique to improve the efficiency of the so-called multiple Bennett acceptance ratio (MBAR) estimator [Shirts and Chodera, J. Chem. Phys.129, 124105 (Year: 2008)] for computing thermodynamic expectations of physical quantities, from samples drawn from Hamiltonian or temperature replica exchange simulations. The methods exploit the Crooks fluctuation theorem[Crooks, J. Stat. Phys.90, 1481 (Year: 1998)] for accurately evaluating the partition functions ratios of neighboring replicas, thus providing an excellent initial guess for the MBAR iterative procedure.

The inhibitors of the Tumor Necrosis Factor-α Converting Enzyme represent promising tools for the treatment of Rheumatoid Arthritis, Multiple Sclerosis and other autoimmune diseases. In this work, using Hamiltonian Replica Exchange Molecular Dynamics simulations and atomistic force field we perform an accurate structural characterization of a group of tartrate-based inhibitors. The simulations highlight a correlation between the conformational landscape in bulk solvent and inhibition potency. Since the structures in bulk solvent are much more compact than the crystallographic bound state, we formulate the hypothesis of a two-step docking mechanism: (i) formation of an intermediate between the compact, hydroxyl exposing conformations in solution and the catalytic zinc ion; (ii) structural rearrangement in the active site of TACE of the zinc-tethered drug in the final binding conformation

Synthetic N-glycosylated CSF114(Glc) and related peptides were proved to be able to recognize specific and high-affinity autoantibodies circulating in blood of relapsing-remitting multiple sclerosis (MS) patients and correlating with disease activity. The effect of these peptides has been linked to the β-turn structure around the minimal epitope Asn(Glc). In this work we performed Hamiltonian replica exchange molecular dynamics simulations on the central heptapeptide fragment of a CSF114(Glc)-derived peptide in water and in a water/hexafluoroacetone mixture, confirming a significant incidence of β-turn structures in both solvents. The structural similarity of the glycosylated and unglycosylated forms in all environments proves that the conformation of the heptapeptide is only marginally affected by the presence of the sugar. Moreover, the presence of a significant amount of bioactive hairpin-like conformations in the water environment suggests a possible use not only in the diagnosis but also in the treatment of MS.

Human VEFGR2 is a receptor tyrosine kinase involved in the angiogenesis process and regarded as an interesting target for the design of anticancer drugs. Its activation/inactivation mechanism is related to conformational changes in its cytoplasmatic kinase domain, involving first among all the αC-helix in N-lobe and the A-loop in C-lobe. Affinity of inhibitors for the active or inactive kinase form could dictate the open or closed conformation of the A-loop, thus making the different conformations of the KDR domain different drug targets in drug discovery. In this view, a detailed knowledge of the conformational landscape of KDR domain is of central relevance to rationalize the efficiency and selectivity of kinase inhibitors. Here, molecular dynamics simulations were used to gain insight into the conformational switching activity of the KDR domain and to identify intermediate conformations between the two limiting active and inactive conformations. Specific energy barriers have been selectively removed to induce, and hence highlight at the atomistic level, the regulation mechanism of the A-loop opening. The proposed strategy allowed to repeatedly observe the escape of the KDR domain from the DFG-out free energy basin and to identify rare intermediate conformations between the DFG-out and the DFG-in structures to be employed in a structure-based drug discovery process.

One of the most important targets in the autoimmune attack in experimental autoimmune encephalomyielitis is the myelin oligodendrocyte glycoprotein (MOG). The complex with demyelinating 8-18C5 antibody was recently resolved by X-ray crystallography, showing a remarkable adhesion of the 101--108 MOG subsequence to the heavy chain of the autoantibody. In this study, we have determined, using replica exchange molecular dynamics methods, the structure of the MOG-derived peptide 101--108 in solution at ambient conditions. According to the simulation, the peptide exhibits, with significant probability, a distorted β-turn structure highly similar to that of the corresponding subsequence in the crystal in complex with 8-18C5 antibody. Such results are found to be fully consistent with circular dichroism spectra of the peptide in solution, suggesting the use of the MOG-derived 101--108 peptide as a potential lead compound for designing decoy targets for the autoimmune attack in multiple sclerosis.

We report on theoretical and experimental data that help elucidate the intimate structural and thermodynamical cooperative mechanisms that are responsible for the inhibitory activity of FK506-related ligands of the FKBP12 protein, a cytosolic enzyme that catalyzes the cis--trans isomerization of prolyl amide bonds and that is involved in immunosuppression and neuronal functioning. Effective FKBP12 ligands are those that mimic in bulk solution the structure of the Tacrolimus natural drug with respect to (i) the orientation of the two carbonyl groups units of the ligand that will form, upon binding, H-bonds with specific residues in the protein binding pocket and (ii) the reduced conformational entropy as in the rigid macrolide Tacrolimus. On this basis, we have rationally designed in silico and synthesized a new compound, Elte421, as a simple variant of an existing FK506-related ligand. The experimental characterization of the synthesized compound via fluorescence quenching shows that Elte421 has a binding affinity for human FKBP12 comparable to that of FK506 natural drug.

The conformational landscape of three FK506-related drugs with disparate inhibition constants is determined in bulk solution using a replica exchange simulation method with solute torsional tempering. Energetic fitness of important drug conformations with respect to the FKBP12 protein is evaluated by molecular mechanics. Results show that the experimental affinity toward peptidyl?prolyl cis?trans isomerase protein (FKBP12) of the analyzed ligands appears to be positively correlated to the observed population of specific chair structures of the drug piperidinic ring in bulk solution. This observation is rationalized on the basis that such structures, stabilized by stereospecific intramolecular hydrophobic interactions, allows the formation of a pair of protein?ligand hydrogen bonds upon binding.

This contribution is aimed at assessing the current knowledge of stacking interactions in proteins focusing on the thermodynamic origin of their strength. First, the energetic nature of p?p interactions is critically discussed by using recent accurate quantum chemical calculations. Then, the intimate relation between energetics and thermodynamics is analyzed via a survey on the last decade computer simulations and experimental data involving interactions of aromatic side chains in solution. The thermodynamics of stacking in protein is further assessed by reviewing studies based on a knowledge-based approach on updated databases of experimental protein structures. Finally the contribution selectively includes some recent authoritative progresses that highlights the paramount importance of stacking in the stabilization and destabilization of protein tertiary structures.

In the framework of single-molecule pulling experiments, the system is typically driven out of equilibrium by a time-dependent external potential V(t) acting on a collective coordinate such that the total Hamiltonian is the sum of V(t) and the Hamiltonian in the absence of external perturbation. Nonequilibrium work theorems such as Jarzynski equality and Crooks fluctuation theorem have been devised to recover free energy differences between states of this extended system. However, one is often more interested in the potential of mean force of the unperturbed Hamiltonian, i.e., in the effective potential dictating the equilibrium distribution of the collective coordinate in the absence of the external potential. In this respect, Hummer and Szabo proposed an algorithm to estimate the desired free energy differences when pulling experiments are performed in only one direction of the process (Proc. Natl. Acad. Sci. USA 2001, 98, 3658). In this paper, we present a potential of mean force estimator of the unperturbed system that exploits the work measured in both forward and backward directions of the process. The method is based on the reweighting technique of Hummer and Szabo and on the Bennett acceptance ratio. Using Brownian-dynamics simulations on a double-well free energy surface, we show that the estimator works satisfactorily in any pulling situation, from nearly equilibrium to strongly dissipative regimes. The method is also applied to the unfolding/refolding process of decaalanine, a system vastly used to illustrate and to test nonequilibrium methodologies. A thorough comparative analysis with another bidirectional potential of mean force estimator (Minh, D. D. L.; Adib, A. B. Phys. Rev. Lett. 2008, 100, 180602) is also presented. The proposed approach is well-suited to recover free energy profiles from single-molecule bidirectional-pulling experiments such as those performed by optical tweezers or atomic force microscopes.

We present non equilibrium molecular dynamics experiments of the unfolding and refolding of a single molecule alanine decapeptide in vacuo subject to a Nosé thermostat. Forward (unfolding) and reverse (refolding) work distribution are numerically calculated for various duration times of the non equilibrium experiments. Crooks theorem is accurately verified for all non equilibrium regimes and the time asymmetry of the process is measured using the recently proposed Jensen-Shannon divergence [E.H. Feng, G. Crooks, Phys. Rev. Lett. 101, 090602 (2008)]. Results on the alanine decapeptide are found similar to recent experimental data on m-RNA molecule in solution, thus evidencing the universal character of the Jensen-Shannon divergence. The patent non Markovianity of the process is rationalized by assuming that the observed forward and reverse distributions can be each described by a combination of two normal distributions satisfying the Crooks theorem, representative of two mutually exclusive linear events. Such bimodal approach reproduces with surprising accuracy the observed non Markovian work distributions. 2010 Elsevier B.V. All rights reserved.

We present the newrelease of the ORAC engine (Procacci et al., J. Comput. Chem. 1997, 18, 1834), a FORTRAN suite to simulate complex biosystems at the atomistic level. The previous release of the ORAC code included multiple time steps integration, smooth particle mesh Ewald method, constant pressure and constant temperature simulations. The present release has been supplemented with the most advanced techniques for enhanced sampling in atomistic systems including replica exchange with solute tempering, metadynamics and steered molecular dynamics. All these computational technologies have been implemented for parallel architectures using the standard MPI communication protocol. ORAC is an open-source program distributed free of charge under the GNU general public license (GPL) at http://www.chim.unifi.it/orac.

In force spectroscopy single-molecule experiments, an individual molecule, usually a polymer, is mechanically stretched by means of an externally controlled driving potential. Typically, the stiffness of this potential is much smaller than the stiffness of the potential of mean force along the molecular extension coordinate. Here we discuss how such a disparity alters the free energy and the reversibility of the driven system, with respect to the pristine molecular system under examination. In particular, by simulating unfolding/refolding experiments of a small protein, we examine the traits of the potential of mean force that are responsible for the dramatic amount of work dissipated in experiments using a soft device. Finally, we show that in bidirectional experiments the free energy of the free molecular system can be easily recovered by appropriate reweighting methods.

n this article we present a potential of mean force estimator based on measurements of the work performed on a system during out of equilibrium realizations of a process. More specifically, the quantities involved in the estimator are the work exponential averages related to the forward and backward directions of the process and the free energy difference between the end states. Such free energy difference can be estimated without resorting to additional methodologies or data, but exploiting the available work measurements in the Bennett acceptance ratio method. Despite the fact that work exponential averages give strongly biased free energy profiles, a simple combination of them, supplied with an accurate estimate of the free energy difference between the end states, provides good free energies, even for fast pulling velocities of the control parameter. Numerical tests have been performed on a deterministic non-Hamiltonian dynamic system (the folding/unfolding process of one alanine deca-peptide) and on a stochastic toy model (a particle which moves into a one-dimensional potential according to Langevin dynamics). In these tests we compare our potential of mean force estimator to the unidirectional Jarzynski equality and to other bidirectional estimators that have appeared in the literature recently.

Multiple Sclerosis (MS) is a highly invalidating autoimmune disease of the Central Nervous System, leading to progressive paralysis and, sometimes, to premature death. One of the potential targets of the autoimmune reaction is the myelin protein MOG (Myelin Oligodendrocyte Glycoprotein). Since the 101?108 fragment of MOG plays a key role in the interaction with the MS-autoantibody 8-18C5, we performed an analysis of the equilibrium conformations of this peptide using the Replica Exchange Molecular Dynamics technique in conjunction with the Generalized Born continuum solvent model. Four variants of the peptide, stabilized by a disulfide bond, were also studied. We found that a significant fraction of the equilibrium population retains the original β-hairpin conformation, and the amount of crystal-like conformations increases in the disulfide-closed analogues. When the equilibrium structures were used in docking simulations with the 8-18C5 autoantibody, we discovered the existence of a docking funnel whose bottom is populated by stable complexes where the peptide occupies the same region of space that was occupied in the crystal. It follows that the MOG 101?108 fragment represents a promising starting point for the design of a drug capable of blocking the 8-18C5 antibody. The molecule may also be used for the development of a diagnostic assay for multiple sclerosis.

We present an approach to the estimate of the potential of mean force along a generic reaction coordinate based on maximum likelihood methods and path-ensemble averages in systems driven far from equilibrium. Following similar arguments, various free energy estimators can be recovered, all providing comparable computational accuracy. The method, applied to the unfolding process of the α-helix form of an alanine decapeptide, gives results in good agreement with thermodynamic integration.

Using a database of 6166 experimental structures taken from the Protein Data Bank, we have studied pair interactions between planar residues (Phe, Tyr, His, Arg, Glu and Asp) in proteins, known as -- interactions. On the basis of appropriate coordinates defining the mutual arrangement of two residues, we have calculated 2-D potentials of mean force aimed at determining the stability of the most probable structures for aromatic--aromatic, aromatic--cation and aromatic--anion bound pairs. Our analysis reveals the thermodynamic relevance and the ubiquity of stacked complexes in proteins.

The generalized Crooks theorem (GCT) for deterministic non-Hamiltonian molecular dynamics simulations [Phys. Rev. E 75, 050101 (2007)] connects the probabilities of nonequilibrium realizations switching the system between two thermodynamic states, to the partition functions of these states. In comparison to the “classical” Crooks nonequilibrium work theorem [J. Stat. Phys. 90, 1481 (1998)], which deals with realizations involving only mechanical work, the GCT also accounts for additional work resulting from changes of the intensive and extensive thermodynamic variables of the system. In this article we present a numerical verification of the GCT using a Lennard-Jones fluid model where two particles are subject to a time-dependent external potential. Moreover, in order to switch the system between different thermodynamic states, the temperature and the pressure (or volume), which are controlled through the Martyna-Tobias-Klein equations of motion [J. Chem. Phys. 101, 4177 (1994)], are also varied externally. The free energy difference between states characterized by different distances of the target particles is evaluated using both a standard methodology (pair radial distribution functions) and the GCT. In order to exploit the various options provided by the GCT approach, i.e., the possibility of temperature/pressure/volume changes during the realizations, the free energy difference is recovered via arbitrary thermodynamic cycles. In all tests, the GCT is quantitatively verified.

We have evaluated the extent to which classical polarizable force fields, based either on the chemical potential equalization principle or on distributed polarizabilities in the framework of the Sum of Interactions Between Fragments Ab initio computed (SIBFA), can reproduce the ab initio polarization energy and the dipole moment of three distinct water oligomers: bifurcated chains, transverse hydrogen-bonded chains, and longitudinal hydrogen-bonded chains of helical shape. To analyze the many-body polarization effect, chains of different size, i.e., from 2 to 12 water monomers, have been considered. Although the dipole moment is a well-defined quantity in both classical polarizable models and quantum mechanical methods, polarization energy can be defined unequivocally only in the former type of approaches. In this study we have used the Kitaura-Morokuma (KM) procedure. Although the KM approach is on the one hand known to overestimate the polarization energy for strongly interacting molecules, on the other hand it can account for the many-body polarization effectively, whereas some other procedures do not. Our data show that, if off-centered lone pair polarizabilities are explicitly represented, classical polarizable force fields can afford a close agreement with the ab initio results, both in terms of polarization energy and in terms of dipole moment.

The Jarzynski identity [C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)] and the Crooks equation [G. E. Crooks, J. Stat. Phys. 90, 1481 (1998)] relate thermodynamic free energy differences to the work done on a system during a collection of nonequilibrium transformations. In the present Rapid Communication we provide generalized versions of these nonequilibrium work theorems, which hold for dissipative transformations where the system may undergo simultaneously mechanical work and pressure-temperature or volume-temperature changes. The proof is valid in the context of dynamic systems that evolve with NPT-based equations of motion according to the Martyna-Tobias-Klein algorithm [Martyna J. Chem. Phys. 101, 4177 (1994)]. An extension of the proof to dynamic systems that evolve through NVT-based equations of motion is also provided. The theorems may be effectively used in non-Hamiltonian molecular dynamics simulations for evaluating Helmholtz or Gibbs free energy differences, or the ratio of partition functions at different temperatures to be eventually used in thermodynamic cycles.

A review of the recent theoretical and computational activity at the Chemistry Department of the University of Firenze in the field of molecular simulations of condensed phases is reported. The topics considered include quantitative methods for accurate free energy calculations, molecular dynamics of liquids and ionic solutions, chemical reactions in solutions, phase transformations and polymerization reactions at high pressures.

The Crooks equation [Eq. (10) in J. Stat. Phys. 90, 1481 (1998)] relates the work done on a system during a nonequilibrium transformation to the free energy difference between the final and the initial state of the transformation. Recently, the authors have derived the Crooks equation for systems in the canonical ensemble thermostatted by the Nose-Hoover or Nose-Hoover chain method [P. Procacci , J. Chem. Phys. 125, 164101 (2006)]. That proof is essentially based on the fluctuation theorem by Evans and Searles [Adv. Phys. 51, 1529 (2002)] and on the equations of motion. Following an analogous approach, the authors derive here the Crooks equation in the context of molecular dynamics simulations of systems in the isothermal-isobaric (NPT) ensemble, whose dynamics is regulated by the Martyna-Tobias-Klein algorithm [J. Chem. Phys. 101, 4177 (1994)]. Their present derivation of the Crooks equation correlates to the demonstration of the Jarzynski identity for NPT systems recently proposed by Cuendet [J. Chem. Phys. 125, 144109 (2006)].

The Crooks equation [Eq. (10) in J. Stat. Phys. 90, 1481 (1998)], originally derived for microscopically reversible Markovian systems, relates the work done on a system during an irreversible transformation to the free energy difference between the final and the initial state of the transformation. In the present work we provide a theoretical proof of the Crooks equation in the context of constant volume, constant temperature steered molecular dynamics simulations of systems thermostated by means of the Nose-Hoover method (and its variant using a chain of thermostats). As a numerical test we use the folding and unfolding processes of decaalanine in vacuo at finite temperature. We show that the distribution of the irreversible work for the folding process is markedly non-Gaussian thereby implying, according to Crooks equation, that also the work distribution of the unfolding process must be inherently non-Gaussian. The clearly asymmetric behavior of the forward and backward irreversible work distributions is a signature of a non-Markovian regime for the folding/unfolding of decaalanine.

Comment on “From subtle to substantial: Role of metal ions on pi-pi interactions”

We propose a new approach for the umbrella sampling method in molecular dynamics simulations of complex systems. An accelerated sampling of the slow degrees of freedom is achieved by generating a single self-adaptive trajectory that tends to span uniformly the reaction coordinate using a time dependent bias potential derived from the preceding history of the system. To show the convergent behavior and the efficiency of the method, we present the free energy surface of alanine dipeptide in water as a function of the backbone dihedral angles.

In the present study we have used molecular dynamics simulations to study the stability of the antiparallel beta-sheet in cellular mouse prion protein (PrPC) and in the D178N mutant. In particular, using the recently developed non-Markovian metadynamics method, we have evaluated the free energy as a function of a reaction coordinate related to the beta-sheet disruption/growth. We found that the antiparallel â-sheet is significantly weaker in the pathogenic D178N mutant than in the wild-type PrPC. The destabilization of PrPC beta-structure in the D178N mutant is correlated to the weakening of the hydrogen bonding network involving the mutated residue, Arg164 and Tyr128 side chains. This in turn indicates that such a network apparently provides a safety mechanism for the unzipping of the antiparallel beta-sheet in the PrPC. We concludethat the antiparallel beta-sheet is likely to undergo disruption rather than growth under pathogenic conditions, in agreement with recent models of the misfolded monomer that assume a parallel alpha-helix.

In the current opinion, the inclusion of polarization response in classical computer simulations is considered as one of the most important and urgent improvements to be implemented in modern empirical potential models. In this work we focus on the capability of polarizable models, based on the pairwise Coulomb interactions, to model systems where strong electric fields enter into play. As shown by Masia, Probst, and Rey (MPR) [in J. Chem. Phys. 121, 7362 (2004)], when a molecule interacts with point charges, polarizable models show underpolarization with respect to ab initio methods. We prove that this underpolarization, clearly related to nonlinear polarization effects, cannot be simply ascribed to the lack of hyperpolarization in the polarizable models, as suggested by MPR. Analysis of the electron-density rearrangement induced on a water molecule by a point charge reveals a twofold level of polarization response. One level involves intramolecular charge transfer on the whole molecular volume, with the related polarization exhibiting a seemingly linear behavior with the external electric field. The other nonlinear polarization level occurs only at strong electric fields and is found to be strictly correlated to the quantum-mechanical nature of the water molecule. The latter type of polarization has a local character, being limited to the space region of the water lone pairs.

Although the cellular monomeric form of the benign prion protein is now well characterized, a model for the monomer of the misfolded conformation (PrPSc) remains elusive. PrPSc quickly aggregates into highly insoluble fibrils making experimental structural characterization very difficult. The tendency to aggregation of PrPSc in aqueous solution implies that the monomer fold must be hydrophobic. Here, by using molecular dynamics simulations, we have studied the cellular mouse prion protein and its D178N pathogenic mutant immersed in a hydrophobic environment (solution of CCl4), to reveal conformational changes and/or local structural weaknesses of the prion protein fold in unfavorable structural and thermodynamic conditions. Simulations in water have been also performed. Although observing in general a rather limited conformation activity in the nanosecond timescale, we have detected a significant weakening of the antiparallel beta-sheet of the D178N mutant in CCl4 and to a less extent in water. No weakening is observed for the native prion protein. The increase of beta-structure in the monomer, recently claimed as evidence for misfolding to PrPSc, has been also observed in this study irrespective of the thermodynamic or structural conditions, showing that this behavior is very likely an intrinsic characteristic of the prion protein fold.

Electronic polarization response in hydrogen-bond clusters and liquid configurations of water and methanol has been studied by density functional theory (DFT) and by a polarizable force field based on the chemical potential equalization (CPE) principle. It has been shown that an accurate CPE parametrization based on isolated molecular properties is not completely transferable to strongly interacting hydrogen-bond clusters with discrepancies between CPE and DFT overall dipole moments as large as 15%. This is due to the lack of intermolecular charge transfer in the standard CPE implementation. A CPE scheme for evaluating the amount of transferred charge has been developed. The charge transfer parameters are determined with the aid of accurate DFT calculations using only hydrogen-bond dimer configurations. The amount of transferred charge is found to be of the order of few hundredths of electrons, as already found in recent studies on hydrogen-bond systems. The parameters of the model are then used, without further adjustment, to different hydrogen-bond clustered forms of water and methanol (oligomer and liquid configurations). In agreement with different approaches proposed in literature for studying charge transfer effects, the transferred charge in hydrogen-bond dimers is found to decrease exponentially with the hydrogen-bond distance. When allowance is made for charge transfer according to the proposed scheme, the CPE dipole moments are found to reproduce satisfactorily the DFT data.

In a recent work [Giese and York J. Chem. Phys. 120, 9903 (2004)] showed that many-body force field models based solely on pairwise Coulomb screening cannot simultaneously reproduce both gas-phase and condensed-phase polarizability limits. In particular, polarizable force fields applied to bifurcated water chains have been demonstrated to be overpolarized with respect to ab initio methods. This behavior was ascribed to the neglect of coupling between many-body exchange and polarization. In the present article we reproduce those results using different ab initio levels of theory and a polarizable model based on the chemical-potential equalization principle. Moreover we show that, when hydrogen-bond (H-bond) forming systems are considered, an additional nonclassical effect, i.e., intermolecular charge transfer, must be taken into account. Such effect, contrarily to that of coupling between many-body exchange and polarization, makes classical polarizable force fields underpolarized. In the case of water at standard conditions, being H-bonded geometries much more probable than the bifurcated ones, intermolecular charge transfer is the dominant effect.

Comment to “Calculation of the dipole moment for polypeptides using the generalized born-electronegativity equalization method: Results in vacuum and continuum-dielectric solvent”

A large set of protein structures resolved by X-ray or NMR techniques has been extracted from the Protein Data Bank and analyzed using statistical methods. In particular, we investigate the interactions between side chains and the interactions between solvent and side chains, pointing out on the possibility of including the solvent as part of a knowledge-based potential. The solvent-residue contacts are accounted for on the basis of the Voronoi's polyhedron analysis. Our investigation confirms the importance of hydrophobic residues in determining the protein stability. We observe that in general hydrophobic-hydrophobic interactions and, more specifically, aromatic-aromatic contacts tend to be increasingly distally separated in the primary sequence of proteins, thus connecting distinct secondary structure elements. A simple relation expressing the dependence of the protein free energy by the number of residues is proposed. Such a relation includes both the residue-residue and the solvent-residue contributions. The former is dominant for large size proteins, whereas for small sizes (number of residues less than 100) the two terms are comparable. Gapless threading experiments show that the solvent-residue knowledge-based potential yields a significant contribution with respect to discriminating the native structure of proteins. Such contribution is important especially for proteins of small size and is similar to that given by the most favorable residue-residue knowledge-based potential referring to hydrophobic-hydrophobic interactions such as isoleucine-leucine. In general, the inclusion of the solvent-residue interaction produces a relevant increase of the free energy gap between the native structures and decoys.

The energetic fitness of histidine in each of its three protonation states has been investigated for NMR-determined protein structures by using molecular mechanics calculations. The protein structures have been taken from the Protein Data Bank (PDB). For the proteins in the database, we generated all isomers, considering all combinations of protonation forms of each histidine. The energy of each isomer has been determined by conjugate gradient minimization using a well-established all-atom force field. We find that, in general, the isomer reported in the PDB is not the most stable isomer. The statistical distribution of isomer energies minus that of the PDB isomer behaves as though the sequence of the histidine forms reported in the PDB was chosen at random. We also show that our molecular mechanics method is a valid approach to predicting the protonation state of histidines buried in the protein core.

This article presents a new ab initio force field for the cofactors of bacterial photosynthesis, namely quinones and bacteriochlorophylls. The parameters has been designed to be suitable for molecular dynamics simulations of photosynthetic proteins by being compatible with the AMBER force field. To our knowledge, this is the first force field for photosynthetic cofactors based on a reliable set of ab initio density functional reference data for methyl bacteriochlorophyll a, methyl bacteriopheophytin a, and of a derivative of ubiquinone. Indeed, the new molecular mechanics force field is able to reproduce very well not only the experimental and ab initio structural properties and the vibrational spectra of the molecules, but also the eigenvectors of the molecular normal modes. For this reason it might also be helpful to understand vibrational spectroscopy results obtained on reaction center proteins.

Structural and dynamical properties of liquid and supercooled liquid m-toluidine are studied by molecular dynamics simulations. Approaching the liquid-glass transition, dynamical heterogeneities, a characteristic common to all supercooled glass formers, are observed. We prove the occurrence of strict correlation between these heterogeneities and the potential energy landscape of the system, expressed in terms of molecule-molecule interactions. A slowing down of the self-diffusive motion of the molecular centers of mass is observed for particular arrangements of pairs of H-bonded molecules. Previous studies on model systems provided evidence of the correlation between dynamical heterogeneities and potential energy landscape, described in terms of the inherent structure of the system. While in this last case the structure is viewed as a collective property of the system, in m-toluidine short-range interactions are sufficient to explain the dynamical behavior in a satisfactory way. This result agrees with the view, supported also by experimental observations, that m-toluidine can be considered as an ensemble of H-bonded subsystems weakly interacting among them.

Comment on “Classical polarizable force fields parametrized from ab initio calculations” [J. Chem. Phys. 117, 1416 (2002)]

A molecular dynamics simulation has been performed to investigate the structure and the dynamics of liquid and supercooled metatoluidine. H-bonding and clustering has been shown to dominate the structure of the liquid. An extensive analysis of the cluster formation and of its relation with the characteristic prepeak of the static structure factor is given. It is shown that molecular association is strongly driven not only by H-bonds, but also by methyl-methyl interactions which favor specific cluster configurations. The liquid-glass transition has been followed through a calculation of the temperature variation of the molar specific heat and a mechanism has been suggested for the structural changes occurring at the phase transition. The librational dynamics of the system has been studied and recent optical Kerr effect measurements have been perfectly reproduced in a wide time regime.

A polarizable electrostatic potential model for classical molecular mechanics is presented. Based on the chemical potential equalization (CPE) principle, the model is developed starting from the original formulation of Mortier, Ghosh, and Shankar [J. Am. Chem. Soc. 108, 4315 (1986)]. Following York and Yang [J. Chem. Phys. 104, 159 (1996)] we present an SP-basis CPE parametrization to describe realistically any sort of molecular system. By fitting ab initio electronic properties, such as dipole moment, polarizability and global molecular hardness of a restricted set of organic molecules, we derive atomic parameters to be applied to a more vast target set of compounds. We show, indeed, that the atomic CPE parameters calculated for the learning set of molecules give reliable values for several electronic properties of various compounds not included in the learning set. The multipole moments obtained by using the proposed CPE parametrization are compared to the results of a fixed charge parametrization like that used by a popular classical molecular mechanics force field, such as AMBER. We show that the fixed charge parametrization can well reproduce only the multipole moments of the molecular conformation or the isomer used for the fit, while being inaccurate when different molecular conformations or isomers are considered. On the contrary, the CPE model realistically reproduces the charge reorganization due to nuclear structural changes of the molecule, such as isomerization or conformational transition. The CPE model has been also tested on various molecular complexes to investigate the polarization response in the case of realistic molecule-molecule interactions. The main result of the paper is the demonstration that the construction of a general polarizable electrostatic force field for classical molecular mechanics is now a viable way.

The nature of intermolecular interactions between aromatic amino acid residues has been investigated by a combination of molecular dynamics and ab initio methods. The potential energy surface of various interacting pairs, including tryptophan, phenilalanine, and tyrosine, was scanned for determining all the relevant local minima by a combined molecular dynamics and conjugate gradient methodology with the AMBER force field. For each of these minima, single-point correlated ab initio calculations of the binding energy were performed. The agreement between empirical force field and ab initio binding energies of the minimum energy structures is excellent. Axomatic-aromatic interactions can be rationalized on the basis of electrostatic and van der Waals interactions, whereas charge transfer or polarization phenomena are small for all intermolecular complexes and, particularly, for stacked structures.

The potential of mean force of interacting aromatic amino acids is calculated using molecular dynamics simulations. The free energy surface is determined in order to study stacking and T-shape competition for phenylalanine-phenylalanine (Phe-Phe), phenylalanine-tyrosine (Phe-Tyr), and tyrosine-tyrosine (Tyr-Tyr) complexes in vacuo, water, carbon tetrachloride, and methanol. Stacked structures are favored in all solvents with the exception of the Tyr-Tyr complex in carbon tetrachloride, where T-shaped structures are also important. The effect of anchoring the two alpha-carbons (C-alpha) at selected distances is investigated. We find that short and large C-alpha-C-alpha distances favor stacked and T-shaped structures, respectively. We analyze a set of 2396 protein structures resolved experimentally. Comparison of theoretical free energies for the complexes to the experimental analogue shows that Tyr-Tyr interaction occurs mainly at the protein surface, while Phe-Tyr and Phe-Phe interactions are more frequent in the hydrophobic protein core. This is confirmed by the Voronoi polyhedron analysis on the database protein structures. As found from the free energy calculation, analysis of the protein database has shown that proximal and distal interacting aromatic residues are predominantly stacked and T-shaped, respectively.

By means of molecular mechanics and ab initio calculations, we show that toluene dimer can assume two different minimum energy structures. Both these arrangements are stacked with the methyl groups being parallel and anti-parallel to each other. Although our findings do not agree with the current opinion that one minimum energy structure is T-shaped, they appear to be consistent with available experiments on jet-cooled toluene.

Modern time resolved spectroscopic techniques, such as off resonant transient birifrangence, or in the frequency domain, such as Fourier transform infrared, constitute a valuable experimental source of information for the dynamical and structural properties of molecular liquids. Many important spectroscopic physical observables can be in principle straightforwardly computed using classical molecular dynamics (MD) simulations. However, MD in its standard implementation has important drawbacks which severely limit its use in the study of light scattering phenomena. In particular, the lack of an adequate representation of the molecular electronic response to the external and local field prevents the accurate simulation of spectroscopic properties. The chemical potential equalization (CPE) method can be used to derive a realistic parameterization for describing electronic polarization effects. The CPE approach is a semi-empirical density functional method. The total molecular energy is expanded to the second order in terms of local atomic electron densities, represented by fluctuating point charges. The empirical parameters of the expansion are atomic electronegativity and hardness which can be fitted onto accurate ab initio calculations of the isolated molecule. CPE is parameterized so as to reproduce global ground state properties of the molecule such as polarizability, ionization potential and electronic affinity as well as local properties such as atomic charges and Fukui indices. From a computational standpoint, the CPE method is simple and inexpensive and can be either straightforwardly implemented in a Car-Parrinello fashion, or used as an analytic tool for computing, a posteriori, spectroscopic properties from trajectories generated by conventional MD simulations, Here we present an extension of the CPE method including atomic dipolar charge distributions and apply it to the simple case of the water monomer and dimer. Very encouraging results are found.

An analysis of structural properties of the tryptophan-histidine (TRP-HIS) pair in vacuo and in various solvents (water, dimethyl sulfoxide, methanol, and carbon tetrachloride) has been done using classical molecular dynamics simulations with atomistic potential models for both solute and solvent molecules. The potential of mean force (PMF) was determined as a function of the TRP centroid-HIS centroid distance and of the angle between the ring planes. We show that T-shaped structures are favored in nonpolar solvents, whereas they are completely destabilized in solvents forming strong hydrogen bonds. Amphiphilic solvents such as methanol and dimethyl sulfoxide destabilize both the T-shaped and stacked arrangements. The results are consistent with previous analysis on protein crystallographic databases, where the stacked structures are found mainly on the protein surface exposed to the solvent, while T-shaped arrangements are preferentially found in the hydrophobic protein core. Finally, comparison of explicit solvent molecular simulation data and continuum solvent model results emphasizes the importance of microsolvation phenomena in shaping the potential of mean force for tryptophan-histidine interactions.

Contrary to the findings of Mülders, Toxvaerd, and Kneller [Phys. Rev. E 58, 6766 (1998)] (MTK), we are unable to discern any difference in the behavior of long chain alkanes simulated by molecular dynamics at constant pressure using either atomic or molecular scaling schemes. This result confirms our previous study [M. Marchi and P. Procacci, J. Chem. Phys. 109, 5194 (1998)] on hydrated proteins published at the same time as the MTK's paper. This Comment indicates that errors in the calculation of the pressure tensor might be responsible for at least a part of the MTKs results.

Erratum: this article contains errata applying to the following content: “Calculation of optical spectra in liquid methanol using molecular dynamics and the chemical potential equalization method”

In this paper we discuss an algorithm for calculation of the temperature-dependent anharmonic correction to the phonon spectrum in atomic and molecular crystals. We show how the equation of motion method can be used to compute corrections of arbitrary perturbation order to the phonon self-energy. Complete analytical expressions up to order λ6 are obtained as an application of the method.

The conformational distribution of glycerol is still an open question both in gas and in liquid phase. Density functional calculations on different conformers of glycerol are reported and compared to the experimental infrared spectra of the gas and of the liquid. The experimental infrared spectra of gas and liquid glycerol are fitted by a linear combination of the single conformer ab initio spectra, obtaining the relative conformer concentrations. For the gas the results are in agreement with electron diffraction experiments and with molecular dynamics simulation data. The conformational distribution of glycerol in liquid phase is less accurate but always indicative. Some results about the role of the intramolecular hydrogen bonding in stabilization and in structural features of the conformers an discussed.

A classical molecular dynamics simulation of liquid benzene is performed, using a potential model which allows for full molecular flexibility. The short range intermolecular radial distribution function is on average reminiscent of the crystalline structure, although practically no preferential orientation can be found for the molecules in the first coordination shell. The average cage lifetime and its vibrational dynamics are obtained from appropriate time correlation functions. The intramolecular vibrations are investigated by calculating the vibrational density of states and the infrared and Raman spectra, achieving an excellent agreement with the experimental data. Finally, the dephasing of the nu (1)(A(1g)) ring breathing mode and of the nu (6)(E-2g) in-plane bending mode is analyzed on the basis of the Kubo dephasing function. For nu (1) mode the Kubo correlation time of 516 fs agrees with the experimental value, and is consistent with a relaxation mechanism involving the cage reorganization. In contrast, nu (6) has a practically pure Lorentzian line shape, with a width of 7.16 cm(-1) in perfect agreement with the experimental value of 7.2 cm(-1).

Empirical force fields for minimum searching in the tryptophan-histidine intermolecular energy surface were used. Fourteen principal minima were identified. For each of these structures the intermolecular energies were computed by using single point correlated ab initio calculation with a split valence and a correlation consistent valence double-zeta basis set. The force field determined complexes have much larger correlated ab initio stabilization energy than those reported in previous studies where a purely ab initio search method was used. The largest stabilization energy was found for a T-shaped complex stabilized by a NH ... N hydrogen bond. Stacked structures with superimposed and parallel-displaced imidazole rings were also found to be very stable.

Combining ab initio and statistical mechanics calculations we have determined the conformational distribution of gas-phase glycerol at different temperatures. The obtained results are consistent with infrared spectroscopy and electron diffraction measurements and are in excellent agreement with previous molecular dynamics simulation data.

Using a model potential function we have performed a molecular dynamics simulation of several static and dynamical properties of glycerol in the crystal, glass and liquid phases. Comparison with available experimental data shows an excellent agreeent and proves the validity of the potential model used. For the calculation of the molar specific heat of the liquid and of the glass we have developed a theoretical approach which takes into account the contributions of the conformational structure energy and of the vibrational energy computed using the Bose-Einstein statistics.

In this paper we compare the polarization response given by two different chemical potential equalization schemes to be applied to molecular dynamics simulations: the standard fluctuating point charge model (FQ) and the atom-atom charge transfer model (AACT). We have tested the transferability of FQ and AACT parameters, fitted to the polarizability of small size alkanes and polyenes, to large size homologues. We show that the FQ scheme is not adequate for the n-alkanes as it strongly overestimates the polarizability tensor components as the number of carbon atoms increases. The FQ approach has been found more predictive for highly conjugated systems like polyenes, although still unsatisfactory. The AACT parameters tuned on ethane are instead perfectly transferable to alkanes of any length and conformation. The AACT scheme satisfactorily reproduces the polarization response also for highly conjugated systems.

We apply the chemical potential equalization (CPE) method to the calculation of the optical spectra in liquid methanol at 298 K and normal pressure. The configurations of the liquid are obtained by conventional molecular dynamics (MD) using a completely flexible all-atoms model. The infrared and Raman spectra are computed a posteriori using a CPE parametrization of methanol calibrated to reproduce the electronic properties of the isolated molecule evaluated with accurate ab initio calculations. The MD/CPE method reproduces correctly the optical spectra in the region of the intermolecular motions. The spectra are discussed and interpreted on the basis of hydrogen bonding structure and dynamics.

An analysis of the conformational properties and hydrogen bonding in the condensed phases of glycerol is reported using the same model as adopted in Part I(Phys. Chem. Chem. Phys., 1999, 1, 871). Structural properties of the liquid and glassy states are analyzed in relation to the molecular backbone conformation of the glycerol molecule. The effects of hydrogen bonding and of temperature on the conformational distribution are analyzed. The structural and dynamical properties of hydrogen bonding in glycerol are also investigated. The results are consistent with available experimental observations and clarify many important and interrelated aspects of the microscopic structure of liquid, glassy and crystalline phases of glycerol.

In this paper, we deal with the handling of the electrostatic forces in complex molecular systems. In particular, we focus on instabilities experienced by reversible multiple time step algorithms when used in conjunction with Ewald summation techniques for periodic systems. We show that energy conservation is negatively affected by the intra-molecular energy term due to electrostatic excluded contacts required by the most frequently used of the modern force fields for biomolecular systems. These effects are due to a non-complete cancellation of the intra-molecular electrostatic energy and forces at intermediate or long time steps.

In this study, we present a new molecular dynamics program for simulation of complex molecular systems. The program, named ORAC, combines state-of-the-art molecular dynamics (MD) algorithms with flexibility in handling different types and sizes of molecules. ORAC is intended for simulations of molecular systems and is specifically designed to treat biomolecules efficiently and effectively in solution or in a crystalline environment. Among its unique features are: (i) implementation of reversible and symplectic multiple time step algorithms (or r-RESPA, reversible reference system propagation algorithm) specifically designed and tuned for biological systems with periodic boundary conditions; (ii) availability for simulations with multiple or single time steps of standard Ewald or smooth particle mesh Ewald (SPME) for computation of electrostatic interactions; and (iii) possibility of simulating molecular systems in a variety of thermodynamic ensembles. We believe that the combination of these algorithms makes ORAC more advanced than other MD programs using standard simulation algorithms.

This paper summarizes, simplifies, and extends some recent developments in the theory of phonon damping in anharmonic crystals. We propose an approach in which the damping is described in terms of two ingredients:m(1) a computed or estimated one-phonon density of states, and (2) average anharmonic couplings between the phonons, fitted to the experimental temperature dependence of the phonon damping. Solid nitrogen is chosen as a test case and the coupling coefficients obtained from the fit are correlated to the three-, four-, and five-phonon couplings computed from a potential model.

Long range electrostatic forces are involved at a fundamental level in many biological phenomena. Their prohibitive computational costs often prevents their correct calculation in molecular dynamics (MD) simulations of biological molecules. In this paper we present a method to handle efficiently and exactly electrostaticinteractions in MD simulations with periodic boundary conditions. Our scheme employs a multiple time step r?RESPA integration algorithm in combination with the Ewald summation technique, and is specifically targeted to simulation of large size complex molecular systems such as solvated proteins. In this approach, the force associated with each particle of the system is partitioned into four components which evolve in time with distinct and increasingly longer time scales. We found that a suitable time scale separation is achieved by subdividing the direct space nonbonded interactions, inclusive of Coulombic and van der Waals contributions, in a short, medium and long range shells. The fastest bonded forces are associated to the component with the shortest timescale, while the reciprocal space nonbonded interaction is included with the medium range direct space forces. Our method is general. It can be straightforwardly implemented for any biomolecular force fields used in condensed phase simulations and can be applied to other complex molecular systems such as molecular liquids and heterogeneous solutions (e.g., micelles, membranes, etc.). We carried out tests on solvated proteins samples containing 7040 and 20627 atoms. Due to the nonlinear scaling of the Ewald computational cost with system size, the performance of our r?RESPA algorithm with respect to single time step algorithms with bond constraints grows with the number of particles. We reach a relative speed?up of 4.2 for a system composed of a solvated photosynthetic reaction center.

In this paper we describe the implementation of a very fast molecular dynamic method to realistically handle electrostatic interactions in simulations of solvated proteins. Our scheme is based on a recently proposed reversible multiple time step (r-RESPA) algorithm and a new modification of the particle mesh Ewald method. While the latter technique provides a fast and accurate representation of the Coulombic interactions for infinite systems, the r-RESPA algorithm exploits a separation of the particle force into components with increasingly longer time scales corresponding to contributions from short-, medium-, and long-range radial shells. By combining the two techniques we are able to reduce considerably the computational cost of molecular dynamics simulation of large biomolecular system without affecting energy conservation and dynamical properties. With respect to single time step simulations employing standard Ewald summation and rigid bond constraints, we obtain a speed-up of about 1 order of magnitude. Finally, our method is about 2.5 times as fast as simulations making use of spherical cutoffs. Since the majority of today biomolecular simulations use spherical cutoffs, we expect that our algorithm will find general applications in the field.

The reversible reference system propagator algorithm (r?RESPA), based on a Trotter factorization of the classical propagator, is tested in a molecular dynamics simulation of solid C60. We show how, with an appropriate subdivision of the interaction potential and with a careful balancing of the integration parameters, one can adopt large time steps and impressively efficient r?RESPA integrators which yield the same dynamics obtained by means of the small time?step Verlet algorithm. The results presented here show that the use of r?RESPA integrators speeds up the simulation by a factor of between 20--40 with respect to the standard Verlet algorithms.

An efficient integration scheme is proposed based on the reversible reference system propagator algorithm (r-RESPA) for the molecular dynamics simulation of molecular solids or liquids at constant pressure. Starting from the Parrinello-Rhaman Lagrangian, the algorithm has been obtained using a Trotter factorization for the classical propagator. This r-RESPA scheme allows us to simulate efficiently the dynamics of completely flexible molecules in the condensed phases, under the constraint of constant pressure. The method may find potential application in the study of mesophases such as nematic or plastic crystals, or in studying the effect of pressure on the spectral properties of solids.

The infrared spectrum of crystalline C60 has been measured in the frequency range 500-20 cm-l at various temperatures between 300 and 10 K. Since no dipole-allowed intramolecular vibration is active below 500 cm-I, the infrared absorption is determined by crystal and quadrupolar interactions. Davydov splittings of several intramolecular vibrations have been observed. The infrared activity of H, modes is discussed. In the lattice region the two translational phonons of Tu symmetry are observed at 41 and 58.5 cm-I at 8 K. Considerations

This paper presents an efficient method for computing high-order corrections to phonon linewidths in crystals. Traditional algorithms involve nested sums over intermediate phonon states, whose computational time grows exponentially with the perturbation order. This nested-sum difficulty is overcome in the present study for the special case of the double-vertex diagrams by writing the corresponding linewidths in a simple form and exploiting a recursive algorithm for weighted phonon densities of states, which requires a time linear in the perturbation order. Using this method, we computed, up to order 10, the linewidths due to the double-vertex diagrams for the α-nitrogen crystal. Our calculation shows that accurate estimates of high-order linewidth corrections can be obtained even for systems as complex as molecular crystals.

A vibrational potential function is proposed for Buckminsterfullerene (C60) in terms of stretching, bending and torsional coordinates. Good agreement of calculated with experimental frequencies is found. On the basis of this potential function a molecular dynamics study of C60 is performed with simulation of the Raman spectrum.

Molecular dynamics simulation studies of rare gases complexes with Buckminsterfullerene are reported in the low-temperature regime. A family of potential parameters is used to describe the C60-rare gas atom interactions and a comparison with different potential models in the case of Xe--C60 interaction is discussed. The effect on the structure and the internal frequencies of C60 is discussed.

The structure of the MCD and absorption spectra of C60 in the region of ? 600 nm is discussed in terms of vibronic activity induced by t1u and h1 modes in the of T1g singlet of lowest energy. Vibronic calculations in the orbital following approach show that t1u vibrations are responsible for the negative while hu for the positive A terms of the MCD spectrum. The assignment of the absorption spectrum is proposed in relation with the MCD results.

In this study is presented a general algorithm for computing Voronoi volumes of atoms of group of atoms in condensed phases. The method is essentially an extension of the Medvedev procedure to allow vertice determination for any Voronoi polyhedron, primitive or with degenerate vertices. The algorithm has been employed for computing time-averaged volumes in the hydrated crystal of met-myoglobin, using the data of a molecular dynamics simulation. The results, compared to previous volume determination in myoglobin, emphasize the fundamental role of solvent structure close to the protein surface in relation to the packing density properties of the residues. Relative volumes fluctuations of myoglobin residues have been found to be correlated to the corresponding mean square displacements from X-ray diffraction studies and to the theoretical hydration energies.

We discuss the double-time Green's-function approach for the phonon-phonon interaction in anharmonic solids. Using a certain interpretation of the phonon operators, the crystal Hamiltonian may be rewritten in a compact form. The equation of motion for the Green's function is tremendously simplified with this formalism. The self-energy can be obtained in a very simple way, to any desired perturbation order. Our results reproduce the corresponding terms of the available literature results, which are limited to the fourth order in the perturbation parameter.

A molecular dynamics study of Xe complexes with fullerenes (C60 and C70) is presented. Two models are used to describe the intermolecular interactions between guest atom and fullerenes, with comparable results. Exohedral complexes may be approximated as van der Waals adducts of Xe with aromatic molecules, with similar equilibrium distances, binding energies and intermolecular frequencies. Endo-complexes are found to be strained due to Xe dimensions and with intermolecular modes much higher in frequency than in the exo-case. The effect of the guest atom on the vibrational dynamics of fullerenes is also investigated.

We have made anharmonic-lattice-dynamics calculations of phonon frequencies and bandwidths up to order λ4 for simple model systems consisting of linear chains of diatomic and triatomic molecules, bound by van der Waals intermolecular forces. For each model we discuss the relevant relaxation mechanisms of the optical modes at the Γ point in terms of the anharmonic phonon-phonon coupling coefficients and of the phonon density of states. Calculated bandwidths and anharmonic shifts are compared to the results of computer simulations for the same models at various temperatures. The comparison shows that bandwidths calculated by the perturbative lattice-dynamics treatment to order λ4 agree with those calculated in the computer simulation, thus showing that the imaginary part of the self-energy expansion converges already to this order. The convergence of the imaginary part of the self-energy is further investigated by a calculation of the contribution of some relevant diagrams of order λ6 and λ8 to the bandwidth of the internal modes. On the contrary, the real part of the self-energy converges much slower and therefore anharmonic shifts calculated by lattice dynamics are different from those obtained by computer simulation.

A vibrational potential function is proposed for C70 in terms of stretching, bending and non-bonded interactions. The calculated frequencies fit the experimental infrared and Raman active vibrations. The force field is discussed and differences and similarities with previous normal mode calculations are pointed out.

The Raman spectrum of a single SO2 crystal in various polarization geometries has been measured at 20 K and discussed in terms of LO and TO components of polar crystal modes. Lattice phonons have been assigned on the basis of their behaviour with polarization. Lattice dynamical calculations on both the internal and lattice phonons are in good agreement with experimental data. Using gas-phase transition dipole moments and a dipole---dipole intermolecular potential the angular dispersion of fundamental vibrons has been calculated. For lattice phonons, a combination of atom---atom and electrostatic potential reproduces not only the phonon frequencies but, for polar phonons, also their LO---TO splitting.

The vibrational relaxation of some internal and external phonons in crystal SO2 as a function of temperature is reported and discussed in terms of depopulation and/or dephasing decay mechanisms. By picosecond CARS spectroscopy the lifetime of the polar ω1 (A1) mode (symmetric stretching) was studied in the temperature range 35--170 K. High-resolution Raman spectroscopy was used to measure the linewidth of the second Davydov component, ω1 (A2), and of lattice modes at 76 and 102 cm ?1 in the same temperature range. The experimental data are analyzed on the basis of current theories on vibrational decay in molecular crystals. Using a recently proposed intermolecular potential for crystal SO2 including atom-atom and electrostatic interactions, the contributions to linewidth from third-order depopulation processes were calculated. In the case of lattice phonons, the calculated bandwidths are in good agreement with experimental results at low temperature. At higher temperature fourth-order decay mechanisms become efficient. As to the internal phonons, calculations show that depopulation has negligible importance for decay. Higher-order contributions to linewidth must be considered in this case. We have estimated the four-phonon intraband scattering and we have found that this mechanism is mostly responsible for the observed behaviour with temperature. It is suggested that fourth-order intraband processes may be important also for the decay of other well isolated phonons in molecular crystals.

Anharmonic lattice dynamics (LD) and molecular dynamics (MD) calculations have been performed for solid CO2 with two different intermolecular potentials. The IR and Raman spectra of the external phonons have been calculated from the MD simulation in the temperature range 25--150 K. The two potential models are shown to be characterized by a different degree of anharmonicity: very broad bands are obtained with the MOSMD potential, while the PRC-1 potential gives narrower spectral lines, in much better agreement with the experimental data. The anharmonic LD calculations agree with the MD results concerning the lineshape; LD instead predicts much larger positive shifts of the phonon frequencies. Such a discrepancy is attributed to the approximation implicit in the perturbative method utilized in LD.

The bandwidths and the anharmonic frequency shifts of the optical lattice phonons of crystalline CO2 are calculated as a function of temperature using three different intermolecular potentials. A previously available potential (MOSMD) which reproduces correctly the structure, energy and harmonic frequencies of the CO2 crystal provides anharmonic, calculated bandwidths about 10 times larger than the experimental values. An improved potential (PRC-1) is obtained by moving inside the O?C bonds the oxygen interaction center and by increasing the molecular quadrupole. This potential predicts bandwidths and shifts of the right order of magnitude but still too large. The third potential (PRC-2) is obtained from the previous one by using a non-collinear distribution of negative charges. The negative charges are localized on exagons around the O?C bond at a distance d = 0.3 Å. This potential reproduces correctly the experimental bandwidths and yields small anharmonic shifts. The temperature dependence of the bandwidths and shifts is correctly reproduced in terms of three-phonon decay processes for the T+g and for the T?g phonons. In the case of the Eg phonon four-phonon processes are required in order to fit the experimental data.