Conformational Energy Research Articles

Simulating protein-ligand binding with neural network potentials.

Drug molecules adopt a range of conformations both in solution and in their protein-bound state. The strain and reduced flexibility of bound drugs can partially counter the intermolecular interactions that drive protein–ligand binding. To make accurate computational predictions of drug binding affinities, computational chemists have attempted to develop efficient empirical models of these interactions, although these methods are not always reliable. Machine learning has allowed the development of highly-accurate neural-network potentials (NNPs), which are capable of predicting the stability of molecular conformations with accuracy comparable to state-of-the-art quantum chemical calculations but at a billionth of the computational cost. Here, we demonstrate that these methods can be used to represent the intramolecular forces of protein-bound drugs within molecular dynamics simulations. These simulations are shown to be capable of predicting the protein–ligand binding pose and conformational component of the absolute Gibbs energy of binding for a set of drug molecules. Notably, the conformational energy for anti-cancer drug erlotinib binding to its target was found to be considerably overestimated by a molecular mechanical model, while the NNP predicts a more moderate value. Although the ANI-1ccX NNP was not trained to describe ionic molecules, reasonable binding poses are predicted for charged ligands, but this method is not suitable for modeling charged ligands in solution.

Overcoming the difficulties of predicting conformational polymorph energetics in molecular crystals via correlated wavefunction methods.

Molecular crystal structure prediction is increasingly being applied to study the solid form landscapes of larger, more flexible pharmaceutical molecules. Despite many successes in crystal structure prediction, van der Waals-inclusive density functional theory (DFT) methods exhibit serious failures predicting the polymorph stabilities for a number of systems exhibiting conformational polymorphism, where changes in intramolecular conformation lead to different intermolecular crystal packings. Here, the stabilities of the conformational polymorphs of o-acetamidobenzamide, ROY, and oxalyl dihydrazide are examined in detail. DFT functionals that have previously been very successful in crystal structure prediction perform poorly in all three systems, due primarily to the poor intramolecular conformational energies, but also due to the intermolecular description in oxalyl dihydrazide. In all three cases, a fragment-based dispersion-corrected second-order Møller-Plesset perturbation theory (MP2D) treatment of the crystals overcomes these difficulties and predicts conformational polymorph stabilities in good agreement with experiment. These results highlight the need for methods which go beyond current-generation DFT functionals to make crystal polymorph stability predictions truly reliable.

Domain-based local pair natural orbital methods within the correlation consistent composite approach.

Ab initio composite approaches have been utilized to model and predict main group thermochemistry within 1 kcal mol-1 , on average, from well-established reliable experiments, primarily for molecules with less than 30 atoms. For molecules of increasing size and complexity, such as biomolecular complexes, composite methodologies have been limited in their application. Therefore, the domain-based local pair natural orbital (DLPNO) methods have been implemented within the correlation consistent composite approach (ccCA) framework, namely DLPNO-ccCA, to reduce the computational cost (disk space, CPU (central processing unit) time, memory) and predict energetic properties such as enthalpies of formation, noncovalent interactions, and conformation energies for organic biomolecular complexes including one of the largest molecules examined via composite strategies, within 1 kcal mol-1 , after calibration with 119 molecules and a set of linear alkanes. © 2019 Wiley Periodicals, Inc.

Fluorescence and Molecular Simulation Studies on the Interaction between Imidazolium-Based Ionic Liquids and Calf Thymus DNA

This work presents a molecular level investigation on the nature and mode of binding between imidazolium-based ionic liquids (ILs) ([Cnbim]Br where n = 2, 4, 6) with calf thymus DNA. This investigation offers valuable insight into the mechanisms of interactions that can affect the structural features of DNA and possibly cause the alteration or inhibition of DNA function. To expedite analysis, the study resorted to using molecular docking and COnductor like Screening MOdel for Real Solvents (COSMO-RS) in conjunction with fluorescence spectroscopic data for confirmation and validation of computational results. Both the fluorescence and docking studies consistently revealed a weak interaction between the two molecules, which corresponded to the binding energy of a stable docking conformation in the range of −5.19 to −7.75 kcal mol−1. As predicted, the rod-like structure of imidazolium-based ILs prefers to bind to the double-helix DNA through a minor groove. Interestingly, the occurrence of T-shape π-π stacking was observed between the amine group in adenine that faces the aromatic ring of imidazole. In addition, data of COSMO-RS for the interaction of individual nucleic acid bases to imidazolium-based ILs affirmed that ILs showed a propensity to bind to different bases, the highest being guanine followed by cytosine, thymine, uracil, and adenine.

β-Arrestin-2 BRET Biosensors Detect Different β-Arrestin-2 Conformations in Interaction with GPCRs.

β-Arrestins are critical regulators of G protein-coupled receptors (GPCRs) that desensitize G protein signaling, promote receptor internalization, and initiate signaling on their own. Recent structural findings indicate that β-arrestins adopt different conformations upon interaction with agonist-activated GPCRs. Here, we established a β-arrestin-2 conformational bioluminescence resonance energy transfer (BRET) sensor composed of the bright Nanoluc BRET donor and the red-shifted CyOFP1 BRET acceptor. The sensor monitors early intramolecular conformational changes of β-arrestin-2 in complex with a wide panel of different class A and class B GPCRs upon agonist activation and with orphan GPCRs known to spontaneously recruit β-arrestin-2. The introduction of the R170E mutant in the β-arrestin-2 sensor allowed the detection of a partially active β-arrestin-2 conformation, which is not dependent on receptor phosphorylation, while the deletion of the β-arrestin-2 finger-loop region detected the "tail-conformation" corresponding to the interaction of β-arrestin with the carboxyl-terminal domain of GPCRs. The new sensors combine the advantages of the BRET technique in terms of sensitivity, robustness, and suitability for real-time measurements with a high responsiveness toward early conformational changes to help to elucidate the different conformational states of β-arrestins associated with GPCR activation in living cells.

Physical Characterization of Inclusion Complexes of Triphenyl Phosphate and Cyclodextrins in Solution.

The goal of this work is to provide physical insights into the formation and stability of inclusion complexes (ICs) in aqueous solution between cyclodextrins (CDs) and a common flame retardant, triphenyl phosphate (TPP). Quantum chemistry calculations reveal the possible energetically favorable geometries of TPP in their 1:1 IC form with α-, β-, and γ-CDs as well as their associated complexation, conformational, and interaction energies. High-resolution mass spectrometry (MS) and tandem MS were used with electrospray ionization to study the soluble ICs formed between TPP and CDs. Successful formation of TPP ICs with both β- and γ-CD in solution was detected in the ratio of 1:1 using high-resolution MS in the positive ion mode. Collision-induced dissociation confirmed the formation of TPP ICs with β- and γ-CDs by generating two product ions, TPP and β- or γ-CD, in both cases. Although quantum chemistry calculations suggest that IC formation with α-CD is energetically possible, an IC with α-CD is not observed in aqueous solution using MS, which aligns with what we also previously observed in the solid state. Since TPP forms stable ICs with β- and γ-CDs both in the solid state and in solution suggests that complexation could be a safer alternative than applying TPP directly to a substrate. In addition, complexation with CDs in solution also opens up new processing methods to create flame-retardant fabrics and foams with TPP.

Perethylated pillar[n]arenes versus pillar[n]arenes: theoretical perspectives.

In this article we obtain the electronic structure, 1H NMR, and the vibrational spectra of perethylated pillar[n]arene (EtPn) (n = 5-8, 10) macrocycles using the ωB97x-based density functional theory. A single cavitand structure endowed with the robust pillar-like architecture was derived for the EtPn (n = 5-7) hosts while for the EtP8 receptor, two conformers one possessing the single cavitand and other double cavitand structure having four hydroquinone monomeric units in each cavity were obtained as the local minima. A double cavitand conformer of the EtP8 was found to be ~ 50kJmol-1 lower in energy than the single cavitand one which was attributed to the increased hydrogen bonding cooperativity and attractive H-H interactions. Further increase of hydroquinone monomers as in EtP10 host engenders only the double cavitand structure, composed of five hydroquinone monomeric units in each cavity. The molecular electrostatic potential investigations along the EtPn series of hosts reveal large electron-rich aromatic moieties compared to unsubstituted pillararenes which suggest the EtPn to be better artificial receptors toward cationic or neutral guests than their unsubstituted analogs. Calculated vibrational frequencies of EtPn (n = 5-8) show the cooperative stretching vibrations in a narrow span of 3241cm-1 to 3228cm-1. Besides, the most deshielded signals from aromatic protons along the EtPn (n = 5-8) series show up near 7.35ppm in their 1H NMR spectra. An increased cooperativity in the lowest energy of the EtP8 conformer further is evident from up-field signal corresponding to methylene protons of ethyl substituent compared with the single cavitand structure. Graphical abstractMolecular electrostatic potential for EtP5 and correlations obtained for δH VsρBCP and ν Vs G(r) plots.

Structural and Dynamic Characterizations Highlight the Deleterious Role of SULT1A1 R213H Polymorphism in Substrate Binding.

Sulfotransferase 1A1 (SULT1A1) is responsible for catalyzing various types of endogenous and exogenous compounds. Accumulating data indicates that the polymorphism rs9282861 (R213H) is responsible for inefficient enzymatic activity and associated with cancer progression. To characterize the detailed functional consequences of this mutation behind the loss-of-function of SULT1A1, the present study deployed molecular dynamics simulation to get insights into changes in the conformation and binding energy. The dynamics scenario of SULT1A1 in both wild and mutated types as well as with and without ligand showed that R213H induced local conformational changes, especially in the substrate-binding loop rather than impairing overall stability of the protein structure. The higher conformational changes were observed in the loop3 (residues, 235–263), turning loop conformation to A-helix and B-bridge, which ultimately disrupted the plasticity of the active site. This alteration reduced the binding site volume and hydrophobicity to decrease the binding affinity of the enzyme to substrates, which was highlighted by the MM-PBSA binding energy analysis. These findings highlight the key insights of structural consequences caused by R213H mutation, which would enrich the understanding regarding the role of SULT1A1 mutation in cancer development and also xenobiotics management to individuals in the different treatment stages.

DFT STUDY OF MOLECULAR STRUCTURE OF 5,10,15,20-TETRAKIS(4’-HALOGENOPHENYL)PORPHYRINS AND THEIR ISOMERS

Conformational manifold of 5,10,15,20-tetrakis(4¢-halogenophenyl)porphyrins (p-4C6H4X-H2P; X=F, Br) was studied by DFT calculations (functional B3LYP). The conformers are different by positions of -4C6H4X groups relative to the macrocyclic core. Relative energies of the conformers of 4C6H4F-H2P were calculated with use of different basis sets in order to find the optimal ratio «quality/ computational cost». According to the results of the calculations, conformers of C2V symmetry are the most energetically favorable. However, relative energies of other conformers are quite low, therefore they should be taken into account at the treatment of the gas-phase electron diffraction (GED) experimental data. The sensitivity of the GED method to structural changes induced by different relative positions of halogenophenyl substituents was examined. Model radial distribution curves f(r) for conformers of p-4C6H4X-H2P, as well as meta- and ortho- isomers (m-4C6H4X-H2P and o-4C6H4X-H2P) were compared. The results of the model studies demonstrate that the bond distances can be reliably determined from experimental data, while the refinement of the positions of halogenophenyl substituents relative to the macrocyclic core is at the limit of the possibilities of the GED method. Meta- and ortho- isomers of 4C6H4X-H2P can be distinguished based on the experimental data, especially in the case of bromine-substituted phenyl groups.

Detecting Functional Dynamics in Proteins with Comparative Perturbed-Ensembles Analysis.

Recent advances have made all-atom molecular dynamics (MD) a powerful tool to sample the conformational energy landscape. There are still however three major challenges in the application of MD to biological systems: accuracy of force field, time scale, and the analysis of simulation trajectories. Significant progress in addressing the first two challenges has been made and extensively reviewed previously. This Account focuses on strategies of analyzing simulation data of biomolecules that also covers ways to properly design simulations and validate simulation results. In particular, we examine an approach named comparative perturbed-ensembles analysis, which we developed to efficiently detect dynamics in protein MD simulations that can be linked to biological functions. In our recent studies, we implemented this approach to understand allosteric regulations in several disease-associated human proteins. The central task of a comparative perturbed-ensembles analysis is to compare two or more conformational ensembles of a system generated by MD simulations under distinct perturbation conditions. Perturbations can be different sequence variations, ligand-binding conditions, and other physical/chemical modifications of the system. Each simulation is long enough (e.g., microsecond-long) to ensure sufficient sampling of the local substate. Then, sophisticated bioinformatic and statistical tools are applied to extract function-related information from the simulation data, including principal component analysis, residue-residue contact analysis, difference contact network analysis (dCNA) based on the graph theory, and statistical analysis of side-chain conformations. Computational findings are further validated with experimental data. By comparing distinct conformational ensembles, functional micro- to millisecond dynamics can be inferred. In contrast, such a time scale is difficult to reach in a single simulation; even when reached for a single condition of a system, it is elusive as to what dynamical motions are related to functions without, for example, comparing free and substrate-bound proteins at the minimum. We illustrate our approach with three examples. First, we discuss using the approach to identify allosteric pathways in cyclophilin A (CypA), a member of a ubiquitous class of peptidyl-prolyl cis-trans isomerase enzymes. By comparing side-chain torsion-angle distributions of CypA in wild-type and mutant forms, we identified three pathways: two are consistent with recent nuclear magnetic resonance experiments, whereas the third is a novel pathway. Second, we show how the approach enables a dynamical-evolution analysis of the human cyclophilin family. In the analysis, both conserved and divergent conformational dynamics across three cyclophilin isoforms (CypA, CypD, and CypE) were summarized. The conserved dynamics led to the discovery of allosteric networks resembling those found in CypA. A residue wise determinant underlying the unique dynamics in CypD was also detected and validated with additional mutational MD simulations. In the third example, we applied the approach to elucidate a peptide sequence-dependent allosteric mechanism in human Pin 1, a phosphorylation-dependent peptidyl-prolyl isomerase. We finally present our outlook of future directions. Especially, we envisage how the approach could help open a new avenue in drug discovery.

Bound Water and Hydroxyproline are the essential contributors to collagen molecular stability: A Computational Analysis

Being the primarily organic phase of bone, collagen type I is an important contributor to bone’s mechanical resistance to fracture. Gaining mechanistic insight into collagen stabilization mechanism is critical to developing new targets to prevent bone fracture. The role of water and hydroxyproline (Hyp) in collagen stability mechanism is still controversial. The aim of this study was to investigate the influences of Hyp and bound water on the collagen molecular stability. Four collagen like-peptide (CLP) models were compared in terms of conformational energies and hydrogen bonding types. CLP1 model represents regular collagen structure without water molecules while CLP2 model represents collagen structure without water and Hyp residue. CLPW1 and CLPW2 are the models of CLP1 and CLP2 with water molecules around them, respectively. Cumulative interpreting of four CLPs models was shed light on the factors influencing collagen stability in the frame of steric energy. Total steric energy was ordered as: CLP2 > CLP1 > CLPW2 > CLPW1, indicating that CLPW1 was the most stable collagen model. On the other hand, CLP2 was the least stable collagen model based on the steric energy comparison. In addition, the hydrogen bonding observed in the four models reveled that water molecules around the models help in binding collagen triple helix through different water bridges since they contributed extra way for binding of triple chains. Moreover, some of the observed water bridges involved directly the presence of Hyp residue. Cumulative results suggested the important role of bound water molecules and Hyp on collagen molecular stability.

Association of functional variants and protein-to-protein physical interactions of human MutY homolog linked with familial adenomatous polyposis and colorectal cancer syndrome.

The human gene MUTYH codes for a DNA glycosylase involved in the repair of oxidative DNA damage. Faulty MUTYH protein activity causes the accumulation of G→T transversions due to unrepaired 8-oxoG:A mismatches. MUTYH germ-line mutations in humans are linked with a recessive form of Familial Adenomatous Polyposis (FAP) and colorectal cancer predisposition. We studied the repair capacity of variants identified in MUTYH-associated polyposis (MAP) patients. MAP is inherited in an autosomal recessive type due to mutations in MUTYH (Y165C, G382D, P54S, A22V, Q63R, G45D, S136P and N43S), indicating that both copies of the gene become inactivated. However, the parents of an individual with an autosomal recessive condition may serve as carriers, each harboring one copy of the mutated gene without showing signs or symptoms of MAP. Six protein partners have been associated with MUTYH, four via direct physical interactions, namely, hMSH6, hPCNA, hRPA1, and hAPEX1. We examined, for the first time, specific interactions of these protein partners with MAP-associated MUTYH mutants using molecular dynamics simulations. The approach provided tools for exploration of the conformational energy landscape accessible to protein partners. The investigation also determined the impact before and after energy minimization of protein-protein interactions and binding affinities of MUTYH wild type and mutant forms, as well as the interactions with other proteins. Taken together, this study provided new insights into the role of MUTYH and its interacting proteins in MAP.

Ramachandran conformational energy maps for disaccharide linkages found in Burkholderia multivorans biofilm polysaccharides

Ramachandran conformational energy maps have been prepared for all of the glycosidic linkages found in the C1576 exopolysaccharide that constitutes the biofilms of the bacterial species Burkholderia multivorans, a member of the Burkholderia cepacian complex that was isolated from a cystic fibrosis patient. This polysaccharide is a rhamnomannan with a tetrasaccharide repeat unit containing two mannose residues and two rhamnose residues, –[3-α-d-Man-(1→2)-α-d-Man-(1→2)-α-d-Rha-(1→3)-α-d-Rha-(1→]n–, where approximately 50% of the rhamnoses are randomly methylated on their O3 hydroxyl groups, further increasing the overall hydrophobicity of the chains. Because of the methylation, the tetrasaccharide repeat unit actually contains six possible linkages. The conformational energy maps are fully adiabatic relaxed maps in which the energy for each (ϕ,ψ) grid point on the map represents the lowest possible energy for the molecule in that conformation, considering all the combinations of the other degrees of freedom, such as hydroxyl orientations. Molecular dynamics simulations were used to verify that these maps indeed describe the conformational dynamics of these linkages. All six linkages were found to be quite restricted in possible ϕ angles, but to exhibit several possible low-energy ψ angles, suggesting that these chains could be quite flexible.

Optimization of synephrine and its vibrational and electronic structures

The optimization of isolated synephrine molecule and its HCl cluster was performed. The potential energy scanning with respect to some dihedral angles at ethylamine side chain was also performed which prefers twelve different low energy structures of single synephrine molecule. The effect of N..H hydrogen bond on geometry and vibrational properties of synephrine in the least energy conformer was studied. The vibrational and electronic structures of all the twelve conformers of synephrine were computed. The experimental FTIR and Raman spectra of synephrine were also recorded in the spectral region 4000–400 cm−1 and 3200-50 cm−1 respectively and correlated with the theoretical spectra of energetically most preferred structure of synephrine. The normal coordinate analysis was employed to scale the overestimated DFT frequencies and to calculate potential energy distributions of vibrtional modes of synephrine. The natural bond orbital analysis was also performed to obtain bond energies and occupancies, HOMO-LUMO and hybridization.

Mechanism of Action of Non-Synonymous Single Nucleotide Variations Associated with α-Carbonic Anhydrase II Deficiency.

Human carbonic anhydrase II (CA-II) is a Zinc (Zn) metalloenzyme responsible for maintenance of acid-base balance within the body through the reversible hydration of CO to produce protons (H) and bicarbonate (BCT). Due to its importance, alterations to the amino acid sequence of the protein as a result of single nucleotide variations (nsSNVs) have detrimental effects on homeostasis. Six pathogenic CA-II nsSNVs, K18E, K18Q, H107Y, P236H, P236R and N252D were identified, and variant protein models calculated using homology modeling. The effect of each nsSNV was analyzed using motif analysis, molecular dynamics (MD) simulations, principal component (PCA) and dynamic residue network (DRN) analysis. Motif analysis identified 11 functionally important motifs in CA-II. RMSD data indicated subtle SNV effects, while PCA analysis revealed that the presence of BCT results in greater conformational sampling and free energy in proteins. DRN analysis showed variant allosteric effects, and the average betweenness centrality (BC) calculations identified Glu117 as the most important residue for communication in CA-II. The presence of BCT was associated with a reduction to Glu117 usage in all variants, suggesting implications for Zn dissociation from the CA-II active site. In addition, reductions to Glu117 usage are associated with increases in the usage of the primary and secondary Zn ligands; His94, His96, His119 and Asn243 highlighting potential compensatory mechanisms to maintain Zn within the active site. Compared to traditional MD simulation investigation, DRN analysis provided greater insights into SNV mechanism of action, indicating its importance for the study of missense mutation effects in proteins and, in broader terms, precision medicine related research.

Time-independent free energies from metadynamics via mean force integration.

Inspired by thermodynamic integration, we propose a method for the calculation of time-independent free energy profiles from history-dependent biased simulations via Mean Force Integration (MFI). MFI circumvents the need for computing the ensemble average of the bias acting on the system c(t) and can be applied to different variants of metadynamics. Moreover, MFI naturally extends to aggregate information obtained from independent metadynamics simulations, allowing to converge free energy surfaces without the need to sample recrossing events in a single continuous trajectory. We validate MFI against one- and two-dimensional analytical potentials and by computing the conformational free energy landscape of ibuprofen in the bulk of its most common crystal phase.

Identification of therapeutic peptide scaffold from tritrpticin family for urinary tract infections using in silico techniques

Antimicrobial peptides (AMPs) like tritrpticins, exhibit non-specific membrane lysis of gram-negative bacteria and can replace antibiotics, combating multi-drug resistance observed in UTI patients. Tritrpticins designated – NT, T1, T2, T3, T5, T7 and T8, were computationally investigated by interaction with Escherichia coli membrane model, mammalian cell toxicity and structural stability to identify a potential drug scaffold for UTI. Initially T3 was eliminated due to low interaction with Escherichia coli membrane model, based on its computed solvation energy. Further, negative support vector machine (SVM) scores revealed non-toxicity of T1, T2, T5, T7 and T8. Finally, at 310 K and varying pH 4.5–9.0, T5 exhibited highest structural stability based on its highest consistency of hydrogen bonds (H-bonds), root mean square deviation (RMSD) and secondary structure profiles along with its lowest conformational free energy. Overall, T5 could be considered a promising peptide drug scaffold to combat UTI.ABBREVIATIONSAMPantimicrobial peptidePBEQPoisson Boltzmann equationH-bondshydrogen bondsMICminimum inhibitory concentrationLD50lethal dose, 50%RMSDroot mean square deviationSVMsupport vector machineUTIurinary tract infectionCommunicated by Ramaswamy H. Sarma

Calculation of Conformational Free Energies with the Focused Confinement Method.

We introduce the focused confinement method, a reaction coordinate-free simulation approach for the calculation of conformational free energies. These are obtained in a series of restrained simulations that transform part of the molecule of interest to independent harmonic oscillators resulting in mixed harmonic-anharmonic states. It is shown that the free energy difference between these mixed states can be readily calculated through the construction of chimeric trajectories. By focusing the confinement to the conformationally active region, the method requires fewer restrained simulations than the traditional confinement method, which eases the treatment of large systems. The accuracy and efficiency of the method is demonstrated for implicitly and explicitly solvated systems.

Ring Puckering Landscapes of Glycosaminoglycan-Related Monosaccharides from Molecular Dynamics Simulations.

The conformational flexibility of the glycosaminoglycans (GAGs) is known to be key in their binding and biological function, for example in regulating coagulation and cell growth. In this work, we employ enhanced sampling molecular dynamics simulations to probe the ring conformations of GAG-related monosaccharides, including a range of acetylated and sulfated GAG residues. We first perform unbiased MD simulations of glucose anomers and the epimers glucuronate and iduronate. These calculations indicate that in some cases, an excess of 15 μs is required for adequate sampling of ring pucker due to the high energy barriers between states. However, by applying our recently developed msesMD simulation method (multidimensional swarm-enhanced sampling molecular dynamics), we were able to quantitatively and rapidly reproduce these ring pucker landscapes. From msesMD simulations, the puckering free energy profiles were then compared for 15 further monosaccharides related to GAGs; this includes to our knowledge the first simulation study of sulfation effects on β-GalNAc ring puckering. For the force field employed, we find that in general the calculated pucker free energy profiles for sulfated sugars were similar to the corresponding unsulfated profiles. This accords with recent experimental studies suggesting that variation in ring pucker of sulfated GAG residues is primarily dictated by interactions with surrounding residues rather than by intrinsic conformational preference. As an exception to this, however, we predict that 4-O-sulfation of β-GalNAc leads to reduced ring rigidity, with a significant lowering in energy of the 1C4 ring conformation; this observation may have implications for understanding the structural basis of the biological function of β-GalNAc-containing glycosaminoglycans such as dermatan sulfate.

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

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