Molecular dynamics simulations of the enzyme catechol-O-methyltransferase: methodological issues

Results from extensive 70 ns all-atom molecular dynamics simulations of catechol-O-methyltransferase (COMT) enzyme are reported. The simulations were performed with explicit TIP3P water and Mg2+ ions. Four different crystal structures of COMT, with and without different ligands, were used. These simulations are among the most extensive of their kind and as such served as a stability test for such simulations. On the methodological side we found that the initial energy minimization procedure may be a crucial step: particular hydrogen bonds may break, and this can initiate an irreversible loss of protein structure that becomes observable in longer time scales of the order of tens of nanoseconds. This has important implications for both molecular dynamics and quantum mechanics-molecular mechanics simulations.     

Molecular dynamics simulations of the enzyme Catechol-O-Methyltransferase: methodological issues1

Results from extensive 70 ns all-atom molecular dynamics simulations of catechol-O-methyltransferase (COMT) enzyme are reported. The simulations were performed with explicit TIP3P water and Mg^2?+ ions. Four different crystal structures of COMT, with and without different ligands, were used. These simulations are among the most extensive of their kind and as such served as a stability test for such simulations. On the methodological side we found that the initial energy minimization procedure may be a crucial step: particular hydrogen bonds may break, and this can initiate an irreversible loss of protein structure that becomes observable in longer time scales of the order of tens of nanoseconds. This has important implications for both molecular dynamics and quantum mechanics–molecular mechanics simulations.     

Computer simulation studies of anisotropic systems. The influence of external fields on the nematic—isotropic transition

The results of Monte Carlo simulations are reported for a model nematogen in the presence of an intense magnetic field. The predictions of a molecular field theory for the system are found to be in reasonable accord with the simulations. It would appear therefore that the critical field for a real nematogen may be unobservably high.     

Comparison of kinetic Monte Carlo and molecular dynamics simulations of diffusion in a model glass former

We present results from kinetic Monte Carlo (KMC) simulations of diffusion in a model glass former. We find that the diffusion constants obtained from KMC simulations have Arrhenius temperature dependence, while the correct behavior, obtained from molecular dynamics simulations, can be super-Arrhenius. We conclude that the discrepancy is due to undersampling of higher-lying local minima in the KMC runs. We suggest that the relevant connectivity of minima on the potential energy surface is proportional to the energy density of the local minima, which determines the "inherent structure entropy." The changing connectivity with potential energy may produce a correlation between dynamics and thermodynamics.     

Path-integral calculations of heavy atom kinetic isotope effects in condensed phase reactions using higher-order Trotter factorizations

A convenient approach to compute kinetic isotope effects (KIEs) in condensed phase chemical reactions is via path integrals (PIs). Usually, the primitive approximation is used in PI simulations, although such quantum simulations are computationally demanding. The efficiency of PI simulations may be greatly improved, if higher-order Trotter factorizations of the density matrix operator are used. In this study, we use a higher-order PI method, in conjunction with mass-perturbation, to compute heavy-atom KIE in the decarboxylation of orotic acid in explicit sulfolane solvent. The results are in good agreement with experiment and show that the mass-perturbation higher-order Trotter factorization provides a practical approach for computing condensed phase heavy-atom KIE.     

Insertion and assembly of membrane proteins via simulation

Interactions of lipids are central to the folding and stability of membrane proteins. Coarse-grained molecular dynamics simulations have been used to reveal the mechanisms of self-assembly of protein/membrane and protein/detergent complexes for representatives of two classes of membrane protein, namely, glycophorin (a simple alpha-helical bundle) and OmpA (a beta-barrel). The accuracy of the coarse-grained simulations is established via comparison with the equivalent atomistic simulations of self-assembly of protein/detergent micelles. The simulation of OmpA/bilayer self-assembly reveals how a folded outer membrane protein can be inserted in a bilayer. The glycophorin/bilayer simulation supports the two-state model of membrane folding, in which transmembrane helix insertion precedes dimer self-assembly within a bilayer. The simulations also suggest that a dynamic equilibrium exists between the glycophorin helix monomer and dimer within a bilayer. The simulated glycophorin helix dimer is remarkably close in structure to that revealed by NMR. Thus, coarse-grained methods may help to define mechanisms of membrane protein (re)folding and will prove suitable for simulation of larger scale dynamic rearrangements of biological membranes.     

Stabilization of a transition-state analogue at the active site of yeast cytosine deaminase: importance of proton transfers

It is believed that the binding of pyrimidin-2-one to cytosine deaminase (CD) leads to the formation of 4-[R]-hydroxyl-3,4-dihydropyrimidine (DHP). Here the formation of transition-state analogue (TSA) at the active site of yeast cytosine deaminase (yCD) is investigated by quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) and free energy simulations. It is shown that DHP may in fact be unstable in the active site and a proton transfer from the Zn hydroxide group to Glu-64 may occur during the nucleophilic attack, leading to an alkoxide-like TSA complex instead. The free energy simulations for the nucleophilic attack process show that the proton transfer from the Zn hydroxide to Glu-64 may play an important role in stabilizing the TSA complex.     

Approximate methods for the fast computation of EPR and ST-EPR spectra. III. Extension of the perturbation approximation to conditions of overmodulatio

The electron paramagnetic resonance (EPR) and saturation transfer electron paramagnetic resonance (ST-EPR) spectra typical of nitroxide spin labels obeying an isotropic, brownian diffusion model are theoretically investigated. The derivative approximation discussed in the first paper of this series (I) is examined and generalized to consider large modulation amplitudes. The approximation considers both molecular and applied modulation effects with the Zeeman modulation frequency acting as a perturbation parameter. The application of the approximation to first and second harmonic simulations is discussed. The observations of I are preserved for first harmonic spectra. The approximation is able to accurately reproduce the in-phase absorption and in-phase-quadrature dispersion signals. In-phase dispersion signals are reproducible at low microwave powers while in-phase-quadrature signals may be reproduced at higher powers. The approximation is of limited utility for second harmonic simulations and is particularly unable to reproduce the in-phase-quadrature absorption signal at this harmonic. The approximation enables a marked reduction in both the core and CPU time required for simulations. Factors of 3 and 4.2 respectively were realized in the simulations of slow isotropic brownian diffusion.     

Structure and Raman spectrum of clavulanic acid in aqueous solution

The calculation of the vibrational Raman spectrum of enzyme-bound beta-lactamase inhibitors may be of help to understand the mechanisms responsible for bacterial drug resistance. Here, we present a study of the solvation structure and the vibrational properties of clavulanate, an important beta-lactamase inhibitor, in aqueous solution as obtained from full quantum and hybrid empirical/quantum molecular dynamics simulations at ambient conditions. The analysis of the vibrational density of states indicates that hybrid empirical/quantum mechanical simulations are able to properly describe the vibrational levels of clavulanate in solution. In addition, we propose a computationally efficient protocol to calculate the vibrational Raman effect for large solute molecules in water, which is able to faithfully reproduce the experimentally recorded clavulanate Raman spectrum and discloses the possibility to employ hybrid simulations to assign the experimental Raman spectra of inhibitors bound to beta-lactamases.     

Molecular dynamics study of taxadiene synthase catalysis

Molecular dynamics (MD) simulations have been performed to study the dynamic behavior of noncovalent enzyme carbocation complexes involved in the cyclization of geranylgeranyl diphosphate to taxadiene catalyzed by taxadiene synthase (TXS). Taxadiene and the observed four side products originate from the deprotonation of carbocation intermediates. The MD simulations of the TXS carbocation complexes provide insights into potential deprotonation mechanisms of such carbocations. The MD results do not support a previous hypothesis that carbocation tumbling is a key factor in the deprotonation of the carbocations by pyrophosphate. Instead water bridges are identified which may allow the formation of side products via multiple proton transfer reactions. A novel reaction path for taxadiene formation is proposed on the basis of the simulations. © 2018 Wiley Periodicals, Inc.     

"All possible steps" approach to the accelerated use of Gillespie's algorithm

Many physical and biological processes are stochastic in nature. Computational models and simulations of such processes are a mathematical and computational challenge. The basic stochastic simulation algorithm was published by Gillespie about three decades ago [J. Phys. Chem. 81, 2340 (1977)]. Since then, intensive work has been done to make the algorithm more efficient in terms of running time. All accelerated versions of the algorithm are aimed at minimizing the running time required to produce a stochastic trajectory in state space. In these simulations, a necessary condition for reliable statistics is averaging over a large number of simulations. In this study the author presents a new accelerating approach which does not alter the stochastic algorithm, but reduces the number of required runs. By analysis of collected data the author demonstrates high precision levels with fewer simulations. Moreover, the suggested approach provides a good estimation of statistical error, which may serve as a tool for determining the number of required runs.     

Molecular dynamics simulation of the LOV2 domain from Adiantum capillus-veneris

The mechanism for signal transduction from the LOV-domains toward the kinase region of phototropin is still not well understood. We have performed molecular dynamics (MD) simulations and CONCOORD calculations on the LOV2 domain of Adiantum capillus-veneris, with the goal to detect possible differences between the two forms of the LOV domain which may not show up in the static crystal structures. Since no such clear differences are found in the MD simulations also, we suggest that the real, biologically active conformation of the LOV domain within the whole phototropin is different from the crystal structure of the isolated LOV domains. The MD simulations do offer, however, insight into details of the dynamics of the dark and illuminated LOV domains, which are discussed in the light of recent experiments.     

Absolute binding free energy calculations of CBClip host–guest systems in the SAMPL5 blind challenge

Herein, we report the absolute binding free energy calculations of CBClip complexes in the SAMPL5 blind challenge. Initial conformations of CBClip complexes were obtained using docking and molecular dynamics simulations. Free energy calculations were performed using thermodynamic integration (TI) with soft-core potentials and Bennett’s acceptance ratio (BAR) method based on a serial insertion scheme. We compared the results obtained with TI simulations with soft-core potentials and Hamiltonian replica exchange simulations with the serial insertion method combined with the BAR method. The results show that the difference between the two methods can be mainly attributed to the van der Waals free energies, suggesting that either the simulations used for TI or the simulations used for BAR, or both are not fully converged and the two sets of simulations may have sampled difference phase space regions. The penalty scores of force field parameters of the 10 guest molecules provided by CHARMM Generalized Force Field can be an indicator of the accuracy of binding free energy calculations. Among our submissions, the combination of docking and TI performed best, which yielded the root mean square deviation of 2.94 kcal/mol and an average unsigned error of 3.41 kcal/mol for the ten guest molecules. These values were best overall among all participants. However, our submissions had little correlation with experiments.     

The self-assembly mechanism of fibril-forming silk-based block copolymers

Triblock copolymers consisting of a silk-based ((Gly-Ala)(3)Gly-Glu) repeat flanked by hydrophilic outer blocks self-assemble into micrometer long fibrils in response to a trigger. Since the exact mechanism of the fibril formation remains unclear, we employ a multiscale modelling approach in combination with rare event simulations to elucidate key processes. Atomistic scale simulations on the silk-based block suggest a mechanism in which a polypeptide prefolded into a β-roll structure docks to the growing end of a fibril through the formation of Glu-Glu sidechain contacts. Subsequently it can slide to the optimal position before water is expelled to form a dry interface between the fibril end and the attaching block copolymer. In addition, we find that the folded state of the silk-based block is further stabilised through interactions with its neighboring block. Templated folding may also play a role in case a partially folded polypeptide attaches. The coarse-grained simulations indicate that the attachment and subsequent sliding is mediated by the hydrophilic flanks in a size dependent manner. The hydrophilic blocks prevent random aggregation and allow growth only at the end of the fibril. Our multiscale approach may be used for other fibril-forming peptides.     

A recent development in computational chemistry: chemical reactions from first principles molecular dynamics simulations

First-principles molecular dynamics simulations have recently been found an effective tool to study a large variety of chemical problems. Finite temperature simulations reveal unique information, including explicit dynamical effects and the evaluation of proper free energy differences. Moreover, dynamics simulations reveal information on the flexibility of molecular systems, and elucidate, often otherwise inaccessible, mechanistic details of chemical reactions. In addition this methodology allows the study of larger, periodic, systems, revealing computationally unique information which may be directly compared to experiments on realistic chemical systems. A variety of examples will be given, although most focus on the important field of catalysis.     

Effect of mutation on enzyme motion in dihydrofolate reductase

Hybrid quantum-classical molecular dynamics simulations of a mutant Escherichia coli dihydrofolate reductase enzyme are presented. Although residue 121 is on the exterior of the enzyme, experimental studies have shown that the mutation of Gly-121 to valine reduces the rate of hydride transfer by a factor of 163. The simulations indicate that the decrease in the hydride transfer rate for the G121V mutant is due to an increase in the free energy barrier. The calculated free energy barrier is higher for the mutant than for the wild-type enzyme by an amount that is consistent with the experimentally observed rate reduction. The calculated transmission coefficients are comparable for the wild-type and mutant enzymes. The simulations suggest that this mutation may interrupt a network of coupled promoting motions proposed to play an important role in DHFR catalysis. This phenomenon has broad implications for protein engineering and drug design.     

Electrostatic potentials in systems periodic in one, two, and three dimensions

We consider the electrostatic potential in a unit cell containing N point charges Q(j) with positions r(j) inside the cell. The cell is replicated periodically in one, two, and three dimensions. The purpose is to give representations for the potential which contain only lattice sums which are absolutely convergent and uniformly convergent in the sampling position r. These representations are derived using variants of the Ewald method and are primarily intended for use in evaluating the accuracy of any algorithm to evaluate electrostatic energies and forces in simulations of dense matter, rather than necessarily for use of themselves in simulations. In reduced dimensionality the Ewald representations can be numerically inefficient and other representations are also provided with careful specification which allows two forms to be used for the potential functions in order to improve numerical performance. These mixed representations may be satisfactory in simulations.     

Anisotropic flow in striped superhydrophobic channels

We report results of dissipative particle dynamics simulations and develop a semi-analytical theory of an anisotropic flow in a parallel-plate channel with two superhydrophobic striped walls. Our approach is valid for any local slip at the gas sectors and an arbitrary distance between the plates, ranging from a thick to a thin channel. It allows us to optimize area fractions, slip lengths, channel thickness, and texture orientation to maximize a transverse flow. Our results may be useful for extracting effective slip tensors from global measurements, such as the permeability of a channel, in experiments or simulations, and may also find applications in passive microfluidic mixing.