A polarizable QM/MM approach to the molecular dynamics of amide groups solvated in water.

The infrared (IR) spectra of polypeptides are dominated by the so-called amide bands. Because they originate from the strongly polar and polarizable amide groups (AGs) making up the backbone, their spectral positions sensitively depend on the local electric fields. Aiming at accurate computations of these IR spectra by molecular dynamics (MD) simulations, which derive atomic forces from a hybrid quantum and molecular mechanics (QM/MM) Hamiltonian, here we consider the effects of solvation in bulk liquid water on the amide bands of the AG model compound N-methyl-acetamide (NMA). As QM approach to NMA we choose grid-based density functional theory (DFT). For the surrounding MM water, we develop, largely based on computations, a polarizable molecular mechanics (PMM) model potential called GP6P, which features six Gaussian electrostatic sources (one induced dipole, five static partial charge distributions) and, therefore, avoids spurious distortions of the DFT electron density in hybrid DFT/PMM simulations. Bulk liquid GP6P is shown to have favorable properties at the thermodynamic conditions of the parameterization and beyond. Lennard-Jones (LJ) parameters of the DFT fragment NMA are optimized by comparing radial distribution functions in the surrounding GP6P liquid with reference data obtained from a "first-principles" DFT-MD simulation. Finally, IR spectra of NMA in GP6P water are calculated from extended DFT/PMM-MD trajectories, in which the NMA is treated by three different DFT functionals (BP, BLYP, B3LYP). Method-specific frequency scaling factors are derived from DFT-MD simulations of isolated NMA. The DFT/PMM-MD simulations with GP6P and with the optimized LJ parameters then excellently predict the effects of aqueous solvation and deuteration observed in the IR spectra of NMA. As a result, the methods required to accurately compute such spectra by DFT/PMM-MD also for larger peptides in aqueous solution are now at hand.

Electrostatics of proteins in dielectric solvent continua. I. An accurate and efficient reaction field description.

We present a reaction field (RF) method which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of an electrostatics calculation with polarizable molecular mechanics (MM) force fields. The method combines an approach originally suggested by Egwolf and Tavan [J. Chem. Phys. 118, 2039 (2003)] with concepts generalizing the Born solution [Z. Phys. 1, 45 (1920)] for a solvated ion. First, we derive an exact representation according to which the sources of the RF potential and energy are inducible atomic anti-polarization densities and atomic shielding charge distributions. Modeling these atomic densities by Gaussians leads to an approximate representation. Here, the strengths of the Gaussian shielding charge distributions are directly given in terms of the static partial charges as defined, e.g., by standard MM force fields for the various atom types, whereas the strengths of the Gaussian anti-polarization densities are calculated by a self-consistency iteration. The atomic volumes are also described by Gaussians. To account for covalently overlapping atoms, their effective volumes are calculated by another self-consistency procedure, which guarantees that the dielectric function ε(r) is close to one everywhere inside the protein. The Gaussian widths σ(i) of the atoms i are parameters of the RF approximation. The remarkable accuracy of the method is demonstrated by comparison with Kirkwood's analytical solution for a spherical protein [J. Chem. Phys. 2, 351 (1934)] and with computationally expensive grid-based numerical solutions for simple model systems in dielectric continua including a di-peptide (Ac-Ala-NHMe) as modeled by a standard MM force field. The latter example shows how weakly the RF conformational free energy landscape depends on the parameters σ(i). A summarizing discussion highlights the achievements of the new theory and of its approximate solution particularly by comparison with so-called generalized Born methods. A follow-up paper describes how the method enables Hamiltonian, efficient, and accurate MM molecular dynamics simulations of proteins in dielectric solvent continua.

Electrostatics of proteins in dielectric solvent continua. II. Hamiltonian reaction field dynamics.

In Paper I of this work [S. Bauer, G. Mathias, and P. Tavan, J. Chem. Phys. 140, 104102 (2014)] we have presented a reaction field (RF) method, which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of polarizable molecular mechanics (MM) force fields. Building upon these results, here we suggest a method for linearly scaling Hamiltonian RF/MM molecular dynamics (MD) simulations, which we call "Hamiltonian dielectric solvent" (HADES). First, we derive analytical expressions for the RF forces acting on the solute atoms. These forces properly account for all those conditions, which have to be self-consistently fulfilled by RF quantities introduced in Paper I. Next we provide details on the implementation, i.e., we show how our RF approach is combined with a fast multipole method and how the self-consistency iterations are accelerated by the use of the so-called direct inversion in the iterative subspace. Finally we demonstrate that the method and its implementation enable Hamiltonian, i.e., energy and momentum conserving HADES-MD, and compare in a sample application on Ac-Ala-NHMe the HADES-MD free energy landscape at 300 K with that obtained in Paper I by scanning of configurations and with one obtained from an explicit solvent simulation.

Piperazine derivatives inhibit PrP/PrP(res) propagation in vitro and in vivo.

Prion diseases are fatal neurodegenerative disorders, which are not curable and no effective treatment exists so far. The major neuropathological change in diseased brains is the conversion of the normal cellular form of the prion protein PrPc(C) into a disease-associated isoform PrP(Sc). PrP(Sc) accumulates into multimeres and fibrillar aggregates, which leads to the formation of amyloid plaques. Increasing evidence indicates a fundamental role of PrP(Sc) species and its aggregation in the pathogenesis of prion diseases, which initiates the pathological cascade and leads to neurodegeneration accompanied by spongiform changes. In search of compounds that have the potential to interfere with PrP(Sc) formation and propagation, we used a cell based assay for the screening of potential aggregation inhibitors. The assay deals with a permanently prion infected cell line that was adapted for a high-throughput screening of a compound library composed of 10,000 compounds (DIVERset 2, ChemBridge). We could detect six different classes of highly potent inhibitors of PrP(Sc) propagation in vitro and identified piperazine derivatives as a new inhibitory lead structure, which increased incubation time of scrapie infected mice.

Spectroscopic polarizable force field for amide groups in polypeptides.

The infrared spectra of polypeptides are dominated by the so-called amide bands. These bands originate from the electrostatically coupled vibrations of the strongly polar amide groups (AGs) making up the polypeptide backbone. Because the AGs are highly polarizable, external electric fields can shift the frequencies of the amide normal modes over wide spectral ranges. The sensitivity to external fields and the strong polarity are the reasons why the shapes of the amide bands can code the structure of the polypeptide backbone. Aiming at a decoding of these band shapes, Schultheis et al. (J. Phys. Chem. B 2008, 112, 12217) have recently suggested a polarizable molecular mechanics (PMM) force field for AGs, which employs field dependent force constants and enables the computation of the amide bands from molecular dynamics simulations. Here we extend and refine this first suggestion of such a PMM force field. The extension rests on the choice of suitable internal coordinates for the AGs and on the inclusion of the complete AG Hessian and of its field dependence. The force field parameters are calculated from density functional theory. The improved quality of the resulting PMM descriptions is demonstrated using very simple examples and an outlook is given.

Highly polar environments catalyze the unfolding of PrP C helix 1.

The first alpha-helix (H1) likely plays an important role in the conversion of the cellular prion protein (PrP(C)) into its pathogenic isoform (PrP(Sc)). In this conversion, H1 may either have to unfold or may represent a site of intermolecular contact. A recent molecular dynamics simulation suggested that H1 can unfold if it is detached from the protein core (Hirschberger et al. in Biophys J 90:3908, 2006). It has been hypothesized that the high dielectric constant epsilon (S) of the bulk water environment facilitates the unfolding of H1. To check this hypothesis, we performed a number of replica exchange molecular dynamics simulations of an H1 peptide in solvents of different epsilon (S). We found that the equilibrium helix fraction in water is less than 40%, in agreement with previous experimental findings, and that the helix unfolds much faster in water than in less polar solvents. The kinetically stabilizing effect of the organic solvents is largely unspecific and correlates well with their dielectric constant epsilon (S).

A polarizable force field for computing the infrared spectra of the polypeptide backbone.

The shapes of the amide bands in the infrared (IR) spectra of proteins and peptides are caused by electrostatically coupled vibrations within the polypeptide backbone and code the structures of these biopolymers. A structural decoding of the amide bands has to resort to simplified models because the huge size of these macromolecules prevents the application of accurate quantum mechanical methods such as density functional theory (DFT). Previous models employed transition-dipole coupling methods that are of limited accuracy. Here we propose a concept for the computation of protein IR spectra, which describes the molecular mechanics (MM) of polypeptide backbones by a polarizable force field of "type II". By extending the concepts of conventional polarizable MM force fields, such a PMM/II approach employs field-dependent parameters not only for the electrostatic signatures of the molecular components but also for the local potentials modeling the stiffness of chemical bonds with respect to elongations, angle deformations, and torsions. Using a PMM/II force field, the IR spectra of the polypeptide backbone can be efficiently calculated from the time dependence of the backbone's dipole moment during a short (e.g., 100 ps) MD simulation by Fourier transformation. PMM/II parameters are derived for harmonic bonding potentials of amide groups in polypeptides from a series of DFT calculations on the model molecule N-methylacetamide (NMA) exposed to homogeneous external electric fields. The amide force constants are shown to vary by as much as 20% for relevant field strengths. As a proof of principle, it is shown that the large solvatochromic effects observed in the IR spectra of NMA upon transfer from the gas phase into aqueous solution are not only excellently reproduced by DFT/MM simulations but are also nicely modeled by the PMM/II approach. The tasks remaining for a proof of practice are specified.

Light-triggered beta-hairpin folding and unfolding.

A light-switchable peptide is transformed with ultrashort pulses from a beta-hairpin to an unfolded hydrophobic cluster and vice versa. The structural changes are monitored by mid-IR probing. Instantaneous normal mode analysis with a Hamiltonian combining density functional theory with molecular mechanics is used to interpret the absorption transients. Illumination of the beta-hairpin state triggers an unfolding reaction that visits several intermediates and reaches the unfolded state within a few nanoseconds. In this unfolding reaction to the equilibrium hydrophobic cluster conformation, the system does not meet significant barriers on the free-energy surface. The reverse folding process takes much longer because it occurs on the time scale of 30 micros. The folded state has a defined structure, and its formation requires an extended search for the correct hydrogen-bond pattern of the beta-strand.

Electrostatics of proteins in dielectric solvent continua. I. Newton's third law marries qE forces.

The authors reformulate and revise an electrostatic theory treating proteins surrounded by dielectric solvent continua [B. Egwolf and P. Tavan, J. Chem. Phys. 118, 2039 (2003)] to make the resulting reaction field (RF) forces compatible with Newton's third law. Such a compatibility is required for their use in molecular dynamics (MD) simulations, in which the proteins are modeled by all-atom molecular mechanics force fields. According to the original theory the RF forces, which are due to the electric field generated by the solvent polarization and act on the partial charges of a protein, i.e., the so-called qE forces, can be quite accurately computed from Gaussian RF dipoles localized at the protein atoms. Using a slightly different approximation scheme also the RF energies of given protein configurations are obtained. However, because the qE forces do not account for the dielectric boundary pressure exerted by the solvent continuum on the protein, they do not obey the principle that actio equals reactio as required by Newton's third law. Therefore, their use in MD simulations is severely hampered. An analysis of the original theory has led the authors now to a reformulation removing the main difficulties. By considering the RF energy, which represents the dominant electrostatic contribution to the free energy of solvation for a given protein configuration, they show that its negative configurational gradient yields mean RF forces obeying the reactio principle. Because the evaluation of these mean forces is computationally much more demanding than that of the qE forces, they derive a suggestion how the qE forces can be modified to obey Newton's third law. Various properties of the thus established theory, particularly issues of accuracy and of computational efficiency, are discussed. A sample application to a MD simulation of a peptide in solution is described in the following paper [M. Stork and P. Tavan, J. Chem. Phys., 126, 165106 (2007).

Identification of anti-prion compounds as efficient inhibitors of polyglutamine protein aggregation in a zebrafish model.

Several neurodegenerative diseases, including Huntington disease (HD), are associated with aberrant folding and aggregation of polyglutamine (polyQ) expansion proteins. Here we established the zebrafish, Danio rerio, as a vertebrate HD model permitting the screening for chemical suppressors of polyQ aggregation and toxicity. Upon expression in zebrafish embryos, polyQ-expanded fragments of huntingtin (htt) accumulated in large SDS-insoluble inclusions, reproducing a key feature of HD pathology. Real time monitoring of inclusion formation in the living zebrafish indicated that inclusions grow by rapid incorporation of soluble htt species. Expression of mutant htt increased the frequency of embryos with abnormal morphology and the occurrence of apoptosis. Strikingly, apoptotic cells were largely devoid of visible aggregates, suggesting that soluble oligomeric precursors may instead be responsible for toxicity. As in nonvertebrate polyQ disease models, the molecular chaperones, Hsp40 and Hsp70, suppressed both polyQ aggregation and toxicity. Using the newly established zebrafish model, two compounds of the N'-benzylidene-benzohydrazide class directed against mammalian prion proved to be potent inhibitors of polyQ aggregation, consistent with a common structural mechanism of aggregation for prion and polyQ disease proteins.

Molecular dynamics simulations indicate a possible role of parallel beta-helices in seeded aggregation of poly-Gln.

The molecular structures of amyloid fibers characterizing neurodegenerative diseases such as Huntington's or transmissible spongiform encephalopathies are unknown. Recently, x-ray diffraction patterns of poly-Gln fibers and electron microscopy images of two-dimensional crystals formed from building blocks of prion rods have suggested that the corresponding amyloid fibers are generated by the aggregation of parallel beta-helices. To explore this intriguing concept, we study the stability of small beta-helices in aqueous solution by molecular dynamics simulations. In particular, for the Huntington aggregation nucleus, which is thought to be formed of poly-Gln polymers, we show that three-coiled beta-helices are unstable at the suggested circular geometries and stable at a triangular shape with 18 residues per coil. Moreover, we demonstrate that individually unstable two-coiled triangular poly-Gln beta-helices become stabilized upon dimerization, suggesting that seeded aggregation of Huntington amyloids requires dimers of at least 36 Gln repeats (or monomers of approximately 54 Gln) for the formation of sufficiently stable aggregation nuclei. An analysis of our results and of sequences occurring in native beta-helices leads us to the proposal of a revised model for the PrP(Sc) aggregation nucleus.

Ultrafast spectroscopy reveals subnanosecond peptide conformational dynamics and validates molecular dynamics simulation.

Femtosecond time-resolved spectroscopy on model peptides with built-in light switches combined with computer simulation of light-triggered motions offers an attractive integrated approach toward the understanding of peptide conformational dynamics. It was applied to monitor the light-induced relaxation dynamics occurring on subnanosecond time scales in a peptide that was backbone-cyclized with an azobenzene derivative as optical switch and spectroscopic probe. The femtosecond spectra permit the clear distinguishing and characterization of the subpicosecond photoisomerization of the chromophore, the subsequent dissipation of vibrational energy, and the subnanosecond conformational relaxation of the peptide. The photochemical cis/trans-isomerization of the chromophore and the resulting peptide relaxations have been simulated with molecular dynamics calculations. The calculated reaction kinetics, as monitored by the energy content of the peptide, were found to match the spectroscopic data. Thus we verify that all-atom molecular dynamics simulations can quantitatively describe the subnanosecond conformational dynamics of peptides, strengthening confidence in corresponding predictions for longer time scales.