Simultaneous parametrization of torsional and third-neighbor interaction terms in force-field development: The LLS-SC algorithm.

The calibration of torsional interaction terms by fitting relative gas-phase conformational energies against their quantum-mechanical values is a common procedure in force-field development. However, much less attention has been paid to the optimization of third-neighbor nonbonded interaction parameters, despite their strong coupling with the torsions. This article introduces an algorithm termed LLS-SC, aimed at simultaneously parametrizing torsional and third-neighbor interaction terms based on relative conformational energies. It relies on a self-consistent (SC) procedure where each iteration involves a linear least-squares (LLS) regression followed by a geometry optimization of the reference structures. As a proof-of-principle, this method is applied to obtain torsional and third-neighbor interaction parameters for aliphatic chains in the context of the GROMOS 53A6 united-atom force field. The optimized parameter set is compared to the original one, which has been fitted manually against thermodynamic properties for small linear alkanes. The LLS-SC implementation is freely available under http://github.com/mssm-labmmol/profiler.

Computer Modeling Explains the Structural Reasons for the Difference in Reactivity of Amine Transaminases Regarding Prochiral Methylketones.

Amine transaminases (ATAs) are pyridoxal-5'-phosphate (PLP)-dependent enzymes that catalyze the transfer of an amino group from an amino donor to an aldehyde and/or ketone. In the past decade, the enzymatic reductive amination of prochiral ketones catalyzed by ATAs has attracted the attention of researchers, and more traditional chemical routes were replaced by enzymatic ones in industrial manufacturing. In the present work, the influence of the presence of an α,ß-unsaturated system in a methylketone model substrate was investigated, using a set of five wild-type ATAs, the (R)-selective from Aspergillus terreus (Atr-TA) and Mycobacterium vanbaalenii (Mva-TA), the (S)-selective from Chromobacterium violaceum (Cvi-TA), Ruegeria pomeroyi (Rpo-TA), V. fluvialis (Vfl-TA) and an engineered variant of V. fluvialis (ATA-256 from Codexis). The high conversion rate (80 to 99%) and optical purity (78 to 99% ee) of both (R)- and (S)-ATAs for the substrate 1-phenyl-3-butanone, using isopropylamine (IPA) as an amino donor, were observed. However, the double bond in the α,ß-position of 4-phenylbut-3-en-2-one dramatically reduced wild-type ATA reactivity, leading to conversions of <10% (without affecting the enantioselectivity). In contrast, the commercially engineered V. fluvialis variant, ATA-256, still enabled an 87% conversion, yielding a corresponding amine with >99% ee. Computational docking simulations showed the differences in orientation and intermolecular interactions in the active sites, providing insights to rationalize the observed experimental results.

On the Effect of the Various Assumptions and Approximations used in Molecular Simulations on the Properties of Bio-Molecular Systems: Overview and Perspective on Issues.

Computer simulations of molecular systems enable structure-energy-function relationships of molecular processes to be described at the sub-atomic, atomic, supra-atomic or supra-molecular level and plays an increasingly important role in chemistry, biology and physics. To interpret the results of such simulations appropriately, the degree of uncertainty and potential errors affecting the calculated properties must be considered. Uncertainty and errors arise from (1) assumptions underlying the molecular model, force field and simulation algorithms, (2) approximations implicit in the interatomic interaction function (force field), or when integrating the equations of motion, (3) the chosen values of the parameters that determine the accuracy of the approximations used, and (4) the nature of the system and the property of interest. In this overview, advantages and shortcomings of assumptions and approximations commonly used when simulating bio-molecular systems are considered. What the developers of bio-molecular force fields and simulation software can do to facilitate and broaden research involving bio-molecular simulations is also discussed.

Drug-Loading Capacity of PAMAM Dendrimers Encapsulating Quercetin Molecules: A Molecular Dynamics Study with the 2016H66 Force Field.

The complexation of quercetin molecules with poly(amidoamine) (PAMAM) dendrimers of generation 0-3 was studied by molecular dynamics simulations. Three main points were addressed: (i) the effect of starting from different initial structures; (ii) the performance of the 2016H66 force field (recently validated in the context of dendrimer simulations) in predicting the experimental drug(quercetin)-loading capacity of PAMAM dendrimers; and (iii) the stability of quercetin-PAMAM complexes and their interactions. Initial structures generated by different restraint protocols led to faster convergence compared to initial structures generated by randomly placing the drug molecules in the simulation box. The simulations yielded meta-stable complexes where the loading numbers have converged to average values and were compared to experimentally obtained values. Once the first meta-stable state was reached, the drug-dendrimer complexes did not deviate significantly throughout the simulation. They were characterized in terms of structural properties, such as the radius of gyration and radial distribution functions. The results suggest that quercetin molecules interact mostly with the internal dendrimer monomers rather than to their surface.

pyPolyBuilder: Automated Preparation of Molecular Topologies and Initial Configurations for Molecular Dynamics Simulations of Arbitrary Supramolecules.

The construction of a molecular topology file is a prerequisite for any classical molecular dynamics simulation. However, the generation of such a file may be very challenging at times, especially for large supramolecules. While many tools are available to provide topologies for large proteins and other biomolecules, the scientific community researching nonbiological systems is not equally well equipped. Here, we present a practical tool to generate topologies for arbitrary supramolecules: The pyPolyBuilder. In addition to linear polymer chains, it also provides the possibility to generate topologies of arbitrary, large, branched molecules, such as, e.g., dendrimers. Furthermore, it also generates reasonable starting structures for simulations of these molecules. pyPolyBuilder is a standalone command-line tool implemented in python. Therefore, it may be easily incorporated in persisting simulation pipelines on any operating systems and with different simulation engines. pyPolyBuilder is freely available on github: https://github.com/mssm-labmmol/pypolybuilder.

gem-Dichlorocyclopropanation of Dicarbonyl Derivatives.

A novel methodology for the 1,1-dichlorocyclopropanation of dicarbonyl conjugated olefins was described. The developed protocol is simple and uses readily accessible starting materials, allowing the isolation of the desired adducts in moderate to excellent yields (up to 99 %). Furthermore, the reaction tolerated scale up to the gram scale; thus highlighting the synthetic potential of this transformation. Control experiments and DFT studies revealed that the reaction proceeded through a Michael-initiated ring-closure process, in which reaction temperature played a crucial role. Finally, these gem-dichlorocyclopropanes were also employed in the preparation of a trisubstituted naphthyl derivative and a diastereoselective reduction was also demonstrated.

Molecular Dynamics Simulations of PAMAM and PPI Dendrimers Using the GROMOS-Compatible 2016H66 Force Field.

A systematic evaluation of the accuracy of the GROMOS-compatible 2016H66 force field in the simulation of dendrimers is performed. More specifically, the poly(amido amine) (PAMAM) and the poly(propyleneimine) (PPI) are considered because of the availability of experimental data and simulation results in the literature. A total of 36 molecular systems are simulated and the radius of gyration, asphericity, density profiles, and the self-diffusion coefficient are monitored in terms of the generation number and pH (low, medium, and high) condition. Overall, the results support the recommendation of the 2016H66 force field for the simulation of dendrimer systems. The natural building-block based strategy adopted in the definition of 2016H66, together with a careful parametrization of the chemical functional groups to reproduce thermodynamic properties in environments of different polarity, and also the ability to accurately reproduce the expected structural and dynamic features of dendrimers, as shown in the present work, make this force field an attractive option for the simulation of such systems.

Solvent effects on the decarboxylation of trichloroacetic acid: insights from ab initio molecular dynamics simulations.

The kinetics of trichloroacetic acid (TCA) decarboxylation strongly depends on the solvent in which it occurs, proceeding faster in polar aprotic solvents compared to protic solvents. In particular, the reaction is known to be fast in DMSO even at room temperature and is rather slow in water even at higher temperatures. In order to understand the role of the solvent in the kinetics of TCA decarboxylation, the present study investigates this reaction using both ab initio molecular dynamics (AIMD) simulations in explicit solvents and static electronic structure calculations with the SMD polarizable continuum model, considering DMSO and water as solvents. Both methodologies yield activation free energies in good agreement with experimental data, however they differ with respect to the reaction profile for the process occurring in water. The simulations suggest that DMSO does not participate chemically in the reaction and that the high reaction rate in DMSO can be explained by differential solvation of the reactant and transition state. In water, a protonation step was observed along the simulation trajectory, indicating chemical participation of the solvent in this case. Moreover, the continuum model has shown to be useful to predict the reaction rates in other solvents, suggesting that reaction rates increase upon decreasing solvent polarity up to the point where the apolar solvents are not able to efficiently screen the strong electrostatic interactions to form the required isolated ionic species.

Insights into CC Chemokine Ligand 2/Chemokine Receptor 2 Molecular Recognition: A Step Forward toward Antichemotactic Agents.

Chemokine ligand 2 (CCL2), also known as monocyte chemoattractant protein 1 (MCP-1), is a chemokine that recruits immune cells to inflammatory sites by interacting with G protein-coupled receptor CCR2. The CCL2/CCR2 axis is also involved in pathological processes such as tumor growth and metastasis and hence is currently considered as an important drug target. CCL2 exists in a dynamic monomer-dimer equilibrium that is modulated by CCR2 binding. We used solution nuclear magnetic resonance (NMR) spectroscopy and molecular dynamics simulations to study the interactions between CCL2 and a sulfopeptide corresponding to the N-terminal sequence of CCR2 (CCR218-31). Peptide binding induced the dissociation of CCL2 into monomers, forming stable CCL2/CCR218-31 complexes. NMR relaxation measurements indicated that residues around the CCR218-31 binding site, which are located at the dimer interface, undergo a complex regime of motions. NMR data were used to construct a three-dimensional structural model of the CCL2/CCR218-31 complex, revealing that CCR218-31 occupies a binding site juxtaposed to the dimer interface, partially replacing monomer-monomer contacts, explaining why CCR218-31 binding weakens the dimer interface and induces dissociation. We found that the main interactions governing receptor binding are highly stable salt bridges with conserved chemokine residues as well as hydrophobic interactions. These data provide new insights into the structure-function relationship of the CCL2-CCR2 interaction and may be helpful for the design of novel antichemotactic agents.

Simulating Bilayers of Nonionic Surfactants with the GROMOS-Compatible 2016H66 Force Field.

Polyoxyethylene glycol alkyl ether amphiphiles (CiEj) are important nonionic surfactants, often used for biophysical and membrane protein studies. In this work, we extensively test the GROMOS-compatible 2016H66 force field in molecular dynamics simulations involving the lamellar phase of a series of CiEj surfactants, namely C12E2, C12E3, C12E4, C12E5, and C14E4. The simulations reproduce qualitatively well the monitored structural properties and their experimental trends along the surfactant series, although some discrepancies remain, in particular in terms of the area per surfactant, the equilibrium phase of C12E5, and the order parameters of C12E3, C12E4, and C12E5. The polar head of the CiEj surfactants is highly hydrated, almost like a single polyethyleneoxide (PEO) molecule at full hydration, resulting in very compact conformations. Within the bilayer, all CiEj surfactants flip-flop spontaneously within tens of nanoseconds. Water-permeation is facilitated, and the bending rigidity is 4 to 5 times lower than that of typical phospholipid bilayers. In line with another recent theoretical study, the simulations show that the lamellar phase of CiEj contains large hydrophilic pores. These pores should be abundant in order to reproduce the comparatively low NMR order parameters. We show that their contour length is directly correlated to the order parameters, and we estimate that they should occupy approximately 7-10% of the total membrane area. Due to their highly dynamic nature (rapid flip-flops, high water permeability, observed pore formation), CiEj surfactant bilayers are found to represent surprisingly challenging systems in terms of modeling. Given this difficulty, the results presented here show that the 2016H66 parameters, optimized independently considering pure-liquid as well as polar and nonpolar solvation properties of small organic molecules, represent a good starting point for simulating these systems.

New insights into flavivirus biology: the influence of pH over interactions between prM and E proteins.

Diseases caused by flaviviruses, such as dengue and zika, are globally recognized as major threats. During infection, a critical point in their replicative cycle is the maturation step, which occurs throughout the cellular exocytic pathway. This step is a pH-dependent process that involves the modification of the viral envelope by converting prM (pre-membrane) into M (membrane) proteins with the release of a "pr peptide". After this reaction, the pr peptides remain bound to the viral envelope while the virions cross the acidic trans-Golgi network, and are released only at neutral pH after secretion of the virus particles. Despite this current knowledge, the molecular basis of the flavivirus maturation step is largely unknown. Here, based on the crystal structure of the dengue pr-E complex ("pr peptide" bound to virus envelope protein) and using molecular dynamics simulations, we found that the pH shift from acidic to neutral yields considerable structural changes in the system. Dynamic cross correlation maps and root mean square deviation analyses revealed that the pr-E junction is clearly unstable under neutral pH. Secondary structure analysis also revealed that the fusion loop region, present in the E protein, is sensitive to pH and tends to unstructure at a neutral environment. Moreover, we found that five residues present in the E protein, Gly102, His244, Thr70, Thr68 and Asn67 are critical to confer stability to the pr-E complex while inside the Golgi apparatus. This work brings details about the dynamical behavior of the pr-E system, helps to better understand the flavivirus biology and may also be of use in the development of novel antiviral strategies.

The mechanism by which P250L mutation impairs flavivirus-NS1 dimerization: an investigation based on molecular dynamics simulations.

The flavivirus non-structural protein 1 (NS1) is a conserved glycoprotein with as yet undefined biological function. This protein dimerizes when inside infected cells or associated to cell membranes but also forms lipid-associated hexamers when secreted to the extracellular space. A single amino acid substitution (P250L) is capable of preventing the dimerization of NS1 resulting in lower virulence and slower virus replication. In this work, based on molecular dynamics simulations of the dengue-2 virus NS1 [Formula: see text]-ladder monomer as a core model, we found that this mutation can induce several conformational changes that importantly affect critical monomer-monomer interactions. Based on additional simulations, we suggest a mechanism by which a highly orchestrated sequence of events propagate the local perturbations around the mutation site towards the dimer interface. The elucidation of such a mechanism could potentially support new strategies for rational production of live-attenuated vaccines and highlights a step forward in the development of novel anti-flavivirus measures.

Origin of the spectral shifts among the early intermediates of the rhodopsin photocycle.

A combined strategy based on the computation of absorption energies, using the ZINDO/S semiempirical method, for a statistically relevant number of thermally sampled configurations extracted from QM/MM trajectories is used to establish a one-to-one correspondence between the structures of the different early intermediates (dark, batho, BSI, lumi) involved in the initial steps of the rhodopsin photoactivation mechanism and their optical spectra. A systematic analysis of the results based on a correlation-based feature selection algorithm shows that the origin of the color shifts among these intermediates can be mainly ascribed to alterations in intrinsic properties of the chromophore structure, which are tuned by several residues located in the protein binding pocket. In addition to the expected electrostatic and dipolar effects caused by the charged residues (Glu113, Glu181) and to strong hydrogen bonding with Glu113, other interactions such as π-stacking with Ala117 and Thr118 backbone atoms, van der Waals contacts with Gly114 and Ala292, and CH/π weak interactions with Tyr268, Ala117, Thr118, and Ser186 side chains are found to make non-negligible contributions to the modulation of the color tuning among the different rhodopsin photointermediates.

Effect of the cosolutes trehalose and methanol on the equilibrium and phase-transition properties of glycerol-monopalmitate lipid bilayers investigated using molecular dynamics simulations.

The influence of the cosolutes trehalose and methanol on the structural, dynamic and thermodynamic properties of a glycerol-1-monopalmitate (GMP) bilayer and on its main transition temperature [Formula: see text] is investigated using atomistic molecular dynamics simulations (600 ns) of a GMP bilayer patch (2 × 8 × 8 lipids) at different temperatures in the range of 302 to 338 K and considering three different cosolute concentrations. Depending on the environment and temperature, these simulations present no or a single GL[Formula: see text]LC, LC[Formula: see text]GL or LC[Formula: see text]ID transition, where LC, GL and ID are the liquid crystal, gel and interdigitated phases, respectively. The trehalose molecules form a coating layer at the bilayer surface, promote the hydrogen-bonded bridging of the lipid headgroups, preserve the interaction of the headgroups with trapped water and induce a slight lateral expansion of the bilayer in the LC phase, observations that may have implications for the phenomenon of anhydrobiosis. However, this cosolute does not affect [Formula: see text] and its dependence on hydration in the concentration range considered. On the other hand, methanol molecules intercalate between the lipid headgroups, promote a lateral expansion of the bilayer in the LC phase and induce a concentration dependent decrease of [Formula: see text], observations that may have implications for the phenomenon of anesthesia. The occurrence of an ID phase in the presence of this cosolute may be viewed as an extreme consequence of lateral expansion. The analysis of the simulations also suggests the existence of two basic conservation principles: (1) the hydrogen-bond saturation principle rests on the observation that for all species present in the different systems, the total numbers of hydrogen-bonds per molecule is essentially constant, the only factor of variability being their distribution among different partners; (2) the densest packing principle rests on the observation that the effective volume per methylene group in the interior of the bilayer is only weakly sensitive to the environment, with values comparable to those for liquid (LC) and solid (ID) alkanes, or intermediate (GL).

gmak: A Parameter-Space Mapping Strategy for Force-Field Calibration.

In the context of classical molecular simulations, the accuracy of a force field is highly influenced by the values of the relevant simulation parameters. In this work, a parameter-space mapping (PSM) workflow is proposed to aid in the calibration of force-field parameters, based mainly on the following features: (i) regular-grid discretization of the search space; (ii) partial sampling of the search-space grid; (iii) training of surrogate models to predict the estimates of the target properties for nonsampled parameter sets; (iv) post hoc interpretation of the results in terms of multiobjective optimization concepts; (v) attenuation of statistical errors achieved via empiric extension of the duration of the simulations; (vi) iterative search-space translation according to a user-defined scalar objective function that measures the accuracy of the force field (e.g., the weighted root-mean-square deviation of the target properties relative to the reference data). This combination of features results in a hybrid of a single- and a multiobjective optimization strategy, allowing for the approximate determination of both a local minimum of the chosen objective function and its neighboring Pareto efficient points. The PSM workflow is implemented in the extensible Python program gmak, which is made available in the Git repository at http://github.com/mssm-labmmol/gmak. Using this implementation, the PSM workflow was tested in a proof-of-concept fashion in the recalibration of the Lennard-Jones parameters of the 3-point Optimal Point Charge (OPC3) water model for compatibility with the GROMOS treatment of nonbonded interactions. The recalibrated model reproduces typical pure-liquid properties with an accuracy similar to the original OPC3 model and represents a significant improvement relative to the Simple Point Charge (SPC) model, which is the official recommendation for simulations using GROMOS force fields.

Reoptimized interaction parameters for the peptide-backbone model compound N-methylacetamide in the GROMOS force field: influence on the folding properties of two beta-peptides in methanol.

Considering N-methylacetamide (NMA) as a model compound, new interaction parameters are developed for the amide function in the GROMOS force field that are compatible with the recently derived 53A6(OXY) parameter set for oxygen-containing chemical functions. The resulting set, referred to as 53A6(OXY+A) , represents an improvement over earlier GROMOS force-field versions in the context of the pure-liquid properties of NMA, including the density, heat of vaporization, dielectric permittivity, self-diffusion constant and viscosity, as well as in terms of the Gibbs hydration free energy of this molecule. Assuming that NMA represents an adequate model compound for the backbone of peptides, 53A6(OXY+A) may be expected to also provide an improved description of polypeptide chains. As an initial test, simulations are reported for two ß-peptides characterized by very different folding properties in methanol. For these systems, earlier force-field versions provided good agreement with experimental NMR data, and the test shows that the improved description achieved in the context of NMA is not accompanied by any deterioration in the representation of the conformational properties of these peptides.

New functionalities in the GROMOS biomolecular simulation software.

Since the most recent description of the functionalities of the GROMOS software for biomolecular simulation in 2005 many new functions have been implemented. In this article, the new functionalities that involve modified forces in a molecular dynamics (MD) simulation are described: the treatment of electronic polarizability, an implicit surface area and internal volume solvation term to calculate interatomic forces, functions for the GROMOS coarse-grained supramolecular force field, a multiplicative switching function for nonbonded interactions, adiabatic decoupling of a number of degrees of freedom with temperature or force scaling to enhance sampling, and nonequilibrium MD to calculate the dielectric permittivity or viscosity. Examples that illustrate the use of these functionalities are given.

Temperature dependence of the dielectric permittivity of acetic acid, propionic acid and their methyl esters: a molecular dynamics simulation study.

For most liquids, the static relative dielectric permittivity is a decreasing function of temperature, because enhanced thermal motion reduces the ability of the molecular dipoles to orient under the effect of an external electric field. Monocarboxylic fatty acids ranging from acetic to octanoic acid represent an exception to this general rule. Close to room temperature, their dielectric permittivity increases slightly with increasing temperature. Herein, the causes for this anomaly are investigated based on molecular dynamics simulations of acetic and propionic acids at different temperatures in the interval 283-363 K, using the GROMOS 53A6(OXY) force field. The corresponding methyl esters are also considered for comparison. The dielectric permittivity is calculated using either the box-dipole fluctuation (BDF) or the external electric field (EEF) methods. The normal and anomalous temperature dependences of the permittivity for the esters and acids, respectively, are reproduced. Furthermore, in the EEF approach, the response of the acids to an applied field of increasing strength is found to present two successive linear regimes before reaching saturation. The low-field permittivity ε, comparable to that obtained using the BDF approach, increases with increasing temperature. The higher-field permittivity ε' is slightly larger, and decreases with increasing temperature. Further analyses of the simulations in terms of radial distribution functions, hydrogen-bonded structures, and diffusion properties suggest that increasing the temperature or the applied field strength both promote a relative population shift from cyclic (mainly dimeric) to extended (chain-like) hydrogen-bonded structures. The lower effective dipole moment associated with the former structures compared to the latter ones provides an explanation for the peculiar dielectric properties of the two acids compared to their methyl esters.

Comparison of the United- and All-Atom Representations of (Halo)alkanes Based on Two Condensed-Phase Force Fields Optimized against the Same Experimental Data Set.

The level of accuracy that can be achieved by a force field is influenced by choices made in the interaction-function representation and in the relevant simulation parameters. These choices, referred to here as functional-form variants (FFVs), include for example the model resolution, the charge-derivation procedure, the van der Waals combination rules, the cutoff distance, and the treatment of the long-range interactions. Ideally, assessing the effect of a given FFV on the intrinsic accuracy of the force-field representation requires that only the specific FFV is changed and that this change is performed at an optimal level of parametrization, a requirement that may prove extremely challenging to achieve in practice. Here, we present a first attempt at such a comparison for one specific FFV, namely the choice of a united-atom (UA) versus an all-atom (AA) resolution in a force field for saturated acyclic (halo)alkanes. Two force-field versions (UA vs AA) are optimized in an automated way using the CombiFF approach against 961 experimental values for the pure-liquid densities ρliq and vaporization enthalpies ΔHvap of 591 compounds. For the AA force field, the torsional and third-neighbor Lennard-Jones parameters are also refined based on quantum-mechanical rotational-energy profiles. The comparison between the UA and AA resolutions is also extended to properties that have not been included as parameterization targets, namely the surface-tension coefficient γ, the isothermal compressibility κT, the isobaric thermal-expansion coefficient αP, the isobaric heat capacity cP, the static relative dielectric permittivity ϵ, the self-diffusion coefficient D, the shear viscosity η, the hydration free energy ΔGwat, and the free energy of solvation ΔGche in cyclohexane. For the target properties ρliq and ΔHvap, the UA and AA resolutions reach very similar levels of accuracy after optimization. For the nine other properties, the AA representation leads to more accurate results in terms of η; comparably accurate results in terms of γ, κT, αP, ϵ, D, and ΔGche; and less accurate results in terms of cP and ΔGwat. This work also represents a first step toward the calibration of a GROMOS-compatible force field at the AA resolution.

Enantiomeric segregation in the gel phase of lipid bilayers.

Enantiospecific interactions within a monoglyceride lipid bilayer are investigated using molecular dynamics simulations. Preferential homochiral interactions are observed in the gel phase, whereas no detectable enantiospecificity is seen in the liquid-crystal phase. On the basis of these results and available experimental data, a mechanism is proposed for the formation of the coagel phase of monoglycerides. Enantiomeric segregation in the gel phase is also discussed in terms of its possible implications for prebiological evolution and membrane raft function.