CHARMM: the biomolecular simulation program

Molecular mechanism by which palmitate inhibits PKR autophosphorylation

PKR (double-stranded RNA-activated protein kinase) is an important component of the innate immunity, antiviral, and apoptotic pathways. Recently, our group found that palmitate, a saturated fatty acid, is involved in apoptosis by reducing the autophosphorylation of PKR at the Thr451 residue; however, the molecular mechanism by which palmitate reduces PKR autophosphorylation is not known. Thus, we investigated how palmitate affects the phosphorylation of the PKR protein at the molecular and biophysical levels. Biochemical and computational studies show that palmitate binds to PKR, near the ATP-binding site, thereby inhibiting its autophosphorylation at Thr451 and Thr446. Mutation studies suggest that Lys296 and Asp432 in the ATP-binding site on the PKR protein are important for palmitate binding. We further confirmed that palmitate also interacts with other kinases, due to the conserved ATP-binding site. A better understanding of how palmitate interacts with the PKR protein, as well as other kinases, could shed light onto possible mechanisms by which palmitate mediates kinase signaling pathways that could have implications on the efficacy of current drug therapies that target kinases.

Stabilization of the predominant disease-causing aldolase variant (A149P) with zwitterionic osmolytes

Hereditary fructose intolerance (HFI) is a disease of carbohydrate metabolism that can result in hyperuricemia, hypoglycemia, liver and kidney failure, coma, and death. Currently, the only treatment for HFI is a strict fructose-free diet. HFI arises from aldolase B deficiency, and the most predominant HFI mutation is an alanine to proline substitution at position 149 (A149P). The resulting aldolase B with the A149P substitution (AP-aldolase) has activity that is <100-fold that of the wild type. The X-ray crystal structure of AP-aldolase at both 4 and 18 °C reveals disordered adjacent loops of the (α/β)(8) fold centered around the substitution, which leads to a dimeric structure as opposed to the wild-type tetramer. The effects of osmolytes were tested for restoration of structure and function. An initial screen of osmolytes (glycerol, sucrose, polyethylene glycol, 2,4-methylpentanediol, glutamic acid, arginine, glycine, proline, betaine, sarcosine, and trimethylamine N-oxide) reveals that glycine, along with similarly structured compounds, betaine and sarcosine, protects AP-aldolase structure and activity from thermal inactivation. The concentration and functional moieties required for thermal protection show a zwitterion requirement. The effects of osmolytes in restoring structure and function of AP-aldolase are described. Testing of zwitterionic osmolytes of increasing size and decreasing fractional polar surface area suggests that osmolyte-mediated AP-aldolase stabilization occurs neither primarily through excluded volume effects nor through transfer free energy effects. These data suggest that AP-aldolase is stabilized by binding to the native structure, and they provide a foundation for developing stabilizing compounds for potential therapeutics for HFI.

The free energy landscape of small molecule unbinding

The spontaneous dissociation of six small ligands from the active site of FKBP (the FK506 binding protein) is investigated by explicit water molecular dynamics simulations and network analysis. The ligands have between four (dimethylsulphoxide) and eleven (5-diethylamino-2-pentanone) non-hydrogen atoms, and an affinity for FKBP ranging from 20 to 0.2 mM. The conformations of the FKBP/ligand complex saved along multiple trajectories (50 runs at 310 K for each ligand) are grouped according to a set of intermolecular distances into nodes of a network, and the direct transitions between them are the links. The network analysis reveals that the bound state consists of several subbasins, i.e., binding modes characterized by distinct intermolecular hydrogen bonds and hydrophobic contacts. The dissociation kinetics show a simple (i.e., single-exponential) time dependence because the unbinding barrier is much higher than the barriers between subbasins in the bound state. The unbinding transition state is made up of heterogeneous positions and orientations of the ligand in the FKBP active site, which correspond to multiple pathways of dissociation. For the six small ligands of FKBP, the weaker the binding affinity the closer to the bound state (along the intermolecular distance) are the transition state structures, which is a new manifestation of Hammond behavior. Experimental approaches to the study of fragment binding to proteins have limitations in temporal and spatial resolution. Our network analysis of the unbinding simulations of small inhibitors from an enzyme paints a clear picture of the free energy landscape (both thermodynamics and kinetics) of ligand unbinding.

Improving the description of salt bridge strength and geometry in a Generalized Born model

The Generalized Born (GB) solvent model is widely used in molecular dynamics simulations because it can be less computationally expensive and it samples conformational changes more efficiently than explicit solvent simulations. Meanwhile, great efforts have been made in the past to improve its precision and accuracy. Previous studies have shown that reducing intrinsic GB radii of some hydrogen atoms would improve AMBER GB-HCT solvent model's accuracy on salt bridges. Here we present our finding that similar correction also shows dramatic improvement for the AMBER GB-OBC solvent model. Potential of mean force and cluster analysis for small peptide replica exchange molecular dynamics simulations suggested that new radii GB simulation with ff99SB/GB-OBC corrected salt bridge strength and achieved significantly higher geometry similarity with TIP3P simulation. Improved performance in 60 ns HIV-1 protease GB simulation further validated this approach for large systems.

RNA polymerase II with open and closed trigger loops: active site dynamics and nucleic acid translocation

RNA polymerase II is the central eukaryotic enzyme in transcription from DNA to RNA. The dynamics of RNA polymerase II is described from molecular-dynamics simulations started from two crystal structures with open and closed trigger loop (TL) forms. Dynamic transitions between neutral and forward translocated states were observed, especially for the downstream DNA duplex. Dynamic rearrangements were also seen in the active site environment, including conformations in which the active site nucleotide assumed a possibly precatalytic conformation in close proximity to the terminal 3'-hydroxyl of the nascent RNA. Because nucleic acid translocation was observed primarily in the simulations with an open TL structure, whereas close approach of the active site nucleotide to the terminal RNA ribose predominantly occurred in the closed TL structure, a modified Brownian ratchet mechanism is proposed whereby thermally driven translocation is only possible with an open TL, and fidelity control and catalysis require TL closing.

Residue-specific side-chain packing determines the backbone dynamics of transmembrane model helices

The transmembrane domains (TMDs) of membrane-fusogenic proteins contain an overabundance of β-branched residues. In a previous effort to systematically study the relation among valine content, fusogenicity, and helix dynamics, we developed model TMDs that we termed LV-peptides. The content and position of valine in LV-peptides determine their fusogenicity and backbone dynamics, as shown experimentally. Here, we analyze their conformational dynamics and the underlying molecular forces using molecular-dynamics simulations. Our study reveals that backbone dynamics is correlated with the efficiency of side-chain to side-chain van der Waals packing between consecutive turns of the helix. Leu side chains rapidly interconvert between two rotameric states, thus favoring contacts to its i±3 and i±4 neighbors. Stereochemical restraints acting on valine side chains in the α-helix force both β-substituents into an orientation where i,i±3 interactions are less favorable than i,i±4 interactions, thus inducing a local packing deficiency at VV3 motifs. We provide a quantitative molecular model to explain the relationship among chain connectivity, side-chain mobility, and backbone flexibility. We expect that this mechanism also defines the backbone flexibility of natural TMDs.

Computing ensembles of transitions from stable states: Dynamic importance sampling

There is an increasing dataset of solved biomolecular structures in more than one conformation and increasing evidence that large-scale conformational change is critical for biomolecular function. In this article, we present our implementation of a dynamic importance sampling (DIMS) algorithm that is directed toward improving our understanding of important intermediate states between experimentally defined starting and ending points. This complements traditional molecular dynamics methods where most of the sampling time is spent in the stable free energy wells defined by these initial and final points. As such, the algorithm creates a candidate set of transitions that provide insights for the much slower and probably most important, functionally relevant degrees of freedom. The method is implemented in the program CHARMM and is tested on six systems of growing size and complexity. These systems, the folding of Protein A and of Protein G, the conformational changes in the calcium sensor S100A6, the glucose-galactose-binding protein, maltodextrin, and lactoferrin, are also compared against other approaches that have been suggested in the literature. The results suggest good sampling on a diverse set of intermediates for all six systems with an ability to control the bias and thus to sample distributions of trajectories for the analysis of intermediate states.

Probing designability via a generalized model of helical bundle geometry

Because the space of folded protein structures is highly degenerate, with recurring secondary and tertiary motifs, methods for representing protein structure in terms of collective physically relevant coordinates are of great interest. By collapsing structural diversity to a handful of parameters, such methods can be used to delineate the space of designable structures (i.e., conformations that can be stabilized with a large number of sequences)-a crucial task for de novo protein design. We first demonstrate this on natural α-helical coiled coils using the Crick parameterization. We show that over 95% of known coiled-coil structures are within  1-Å C(α) root mean square deviation of a Crick-ideal backbone. Derived parameters show that natural geometric space of coiled coils is highly restricted and can be represented by "allowed" conformations amidst a potential continuum of conformers. Allowed structures have (1) restricted axial offsets between helices, which differ starkly between parallel and anti-parallel structures; (2) preferred superhelical radii, which depend linearly on the oligomerization state; (3) pronounced radius-dependent a- and d-position amino acid propensities; and (4) discrete angles of rotation of helices about their axes, which are surprisingly independent of oligomerization state or orientation. In all, we estimate the space of designable coiled-coil structures to be reduced at least 160-fold relative to the space of geometrically feasible structures. To extend the benefits of structural parameterization to other systems, we developed a general mathematical framework for parameterizing arbitrary helical structures, which reduces to the Crick parameterization as a special case. The method is successfully validated on a set of non-coiled-coil helical bundles, frequent in channels and transporter proteins, which show significant helix bending but not supercoiling. Programs for coiled-coil parameter fitting and structure generation are provided via a web interface at http://www.gevorggrigoryan.com/cccp/, and code for generalized helical parameterization is available upon request.

Conformational properties of methyl β-maltoside and methyl α- and β-cellobioside disaccharides

An investigation of the conformational properties of methyl β-maltoside, methyl α-cellobioside, and methyl β-cellobioside disaccharides using NMR spectroscopy and molecular dynamics (MD) techniques, is presented. Emphasis is placed on validation of a recently presented force field for hexopyranose disaccharides followed by elucidation of the conformational properties of two different types of glycosidic linkages, α-(1 → 4) and β-(1 → 4). Both gas-phase and aqueous-phase simulations are performed to gain insight into the effect of solvent on the conformational properties. A number of transglycosidic J-coupling constants and proton-proton distances are calculated from the simulations and are used to identify the percent sampling of the three glycosidic conformations (syn, anti-φ, and anti-ψ) and, in turn, describe the flexibility around the glycosidic linkage. The results show the force field to be in overall good agreement with experiment, although some very small limitations are evident. Subsequently, a thorough hydrogen bonding analysis is performed to obtain insights into the conformational properties of the disaccharides. In methyl β-maltoside, competition between HO2'-O3 intramolecular hydrogen bonding and intermolecular hydrogen bonding of those groups with solvent leads to increased sampling of syn, anti-φ, and anti-ψ conformations and better agreement with NMR J-coupling constants. In methyl α- and β-cellobioside, O5'-HO6 and HO2'-O3 hydrogen bonding interactions are in competition with intermolecular hydrogen bonding involving the solvent molecules. This competition leads to retention of the O5'-HO3 hydrogen bond and increased sampling of the syn region of the φ/ψ map. Moreover, glycosidic torsions are correlated to the intramolecular hydrogen bonding occurring in the molecules. The present results verify that in the β-(1 → 4)-linkage intramolecular hydrogen bonding in the aqueous phase is due to the decreased ability of water to successfully compete for the O5' and HO3 hydrogen bonding moieties, in contrast to that occurring between the O5' and HO6 atoms in this α-(1 → 4)-linkage.

Protegrin-1 orientation and physicochemical properties in membrane bilayers studied by potential of mean force calculations

Protegrin-1 (PG-1) belongs to the family of antimicrobial peptides. It interacts specifically with the membrane of a pathogen and kills the pathogen by releasing its cellular contents. To fully understand the energetics governing the orientation of PG-1 in different membrane environments and its effects on the physicochemical properties of the peptide and membrane bilayers, we have performed the potential of mean force (PMF) calculations as a function of its tilt angle at four distinct rotation angles in explicit membranes composed of either DLPC (1,2-dilauroylphosphatidylcholine) or POPC (1-palmitoyl-2-oleoylphosphatidylcholine) lipid molecules. The resulting PMFs in explicit lipid bilayers were then used to search for the optimal hydrophobic thickness of the EEF1/IMM1 implicit membrane model in which a two-dimensional PMF in the tilt and rotation space was calculated. The PMFs in explicit membrane systems clearly reveal that the energetically favorable tilt angle is affected by both the membrane hydrophobic thickness and the PG-1 rotation angle. Local thinning of the membrane around PG-1 is observed upon PG-1 tilting. The thinning is caused by both hydrophobic mismatch and arginine-lipid head group interactions. The two-dimensional PMF in the implicit membrane is in good accordance with those from the explicit membrane simulations. The ensemble-averaged Val16 (15)N and (13)CO chemical shifts weighted by the two-dimensional PMF agree fairly well with the experimental values, suggesting the importance of peptide dynamics in calculating such ensemble properties for direct comparison with experimental observables.

Development of CHARMM polarizable force field for nucleic acid bases based on the classical Drude oscillator model

A polarizable force field for nucleic acid bases based on the classical Drude oscillator model is presented. Parameter optimization was performed to reproduce crystallographic geometries, crystal unit cell parameters, heats of sublimation, vibrational frequencies and assignments, dipole moments, molecular polarizabilities and quantum mechanical base-base and base-water interaction energies. The training and validation data included crystals of unsubstituted and alkyl-substituted adenine, guanine, cytosine, uracil, and thymine bases, hydrated crystals, and hydrogen bonded base pairs. Across all compounds, the RMSD in the calculated heats of sublimation is 4.1%. This equates to an improvement of more than 2.5 kcal/mol in accuracy compared to the nonpolarizable CHARMM27 force field. However, the level of agreement with experimental molecular volume decreased from 1.7% to 2.1% upon moving from the nonpolarizable to the polarizable model. The representation of dipole moments is significantly improved with the Drude polarizable force field. Unlike in additive force fields, there is no requirement for the gas-phase dipole moments to be overestimated, illustrating the ability of the Drude polarizable force field to treat accurately differently dielectric environments and indicating the improvements in the electrostatic model. Validation of the model was performed on the basis of the calculation of the gas-phase binding enthalpies of base pairs obtained via potential of mean force calculations; the additive and polarizable models both performed satisfactorily with average differences of 0.2 and 0.9 kcal/mol, respectively, and rms differences of 1.3 and 1.7 kcal/mol, respectively. Overall, considering the number of significant improvements versus the additive CHARMM force field, the incorporation of explicit polarizability into the force field for nucleic acid bases represents an additional step toward accurate computational modeling of biological systems.

CHARMM additive all-atom force field for glycosidic linkages in carbohydrates involving furanoses

Presented is an extension of the CHARMM additive carbohydrate all-atom force field to enable modeling of polysaccharides containing furanose sugars. The new force field parameters encompass 1 ↔ 2, 1 → 3, 1 → 4, and 1 → 6 pyranose-furanose linkages and 2 → 1 and 2 → 6 furanose-furanose linkages, building on existing hexopyranose and furanose monosaccharide parameters. The model compounds were chosen to be monomers or glycosidic-linked dimers of tetrahydropyran (THP) and tetrahydrofuran (THF) as to contain the key atoms in full carbohydrates. Target data for optimization included two-dimensional quantum mechanical (QM) potential energy scans of the Φ/Ψ glycosidic dihedral angles, with geometry optimization at the MP2/6-31G(d) level followed by MP2/cc-pVTZ single-point energies. All possible chiralities of the model compounds at the linkage carbons were considered, and for each geometry, the THF ring was constrained to the favorable south or north conformations. Target data also included QM vibrational frequencies and pair interaction energies and distances with water molecules. Force field validation included comparison of computed crystal properties, aqueous solution densities, and NMR J-coupling constants to experimental reference values. Simulations of infinite crystals showed good agreement with experimental values for intramolecular geometries as well as for crystal unit cell parameters. Additionally, aqueous solution densities and available NMR data were reproduced to a high degree of accuracy, thus validating the hierarchically optimized parameters in both crystalline and aqueous condensed phases. The newly developed parameters allow for the modeling of linear, branched, and cyclic pyranose/furanose polysaccharides both alone and in heterogeneous systems including proteins, nucleic acids, and/or lipids when combined with existing additive CHARMM biomolecular force fields.

Unfolding and translocation pathway of substrate protein controlled by structure in repetitive allosteric cycles of the ClpY ATPase

Clp ATPases are ring-shaped AAA+ motors in the degradation pathway that perform critical actions of unfolding and translocating substrate proteins (SPs) through narrow pores to deliver them to peptidase components. These actions are effected by conserved diaphragm-forming loops found in the central channel of the Clp ATPase hexamer. Conformational changes, that take place in the course of repetitive ATP-driven cycles, result in mechanical forces applied by the central channel loops onto the SP. We use coarse-grained simulations to elucidate allostery-driven mechanisms of unfolding and translocation of a tagged four-helix bundle protein by the ClpY ATPase. Unfolding is initiated at the tagged C-terminal region via an obligatory intermediate. The resulting nonnative conformation is competent for translocation, which proceeds on a different time scale than unfolding and involves sharp stepped transitions. Completion of the translocation process requires assistance from the ClpQ peptidase. These mechanisms contrast nonallosteric mechanical unfolding of the SP. In atomic force microscopy experiments, multiple unfolding pathways are available and large mechanical forces are required to unravel the SP relative to those exerted by the central channel loops of ClpY. SP threading through a nonallosteric ClpY nanopore involves simultaneous unfolding and translocation effected by strong pulling forces.

RNA and protein 3D structure modeling: similarities and differences

In analogy to proteins, the function of RNA depends on its structure and dynamics, which are encoded in the linear sequence. While there are numerous methods for computational prediction of protein 3D structure from sequence, there have been very few such methods for RNA. This review discusses template-based and template-free approaches for macromolecular structure prediction, with special emphasis on comparison between the already tried-and-tested methods for protein structure modeling and the very recently developed “protein-like” modeling methods for RNA. We highlight analogies between many successful methods for modeling of these two types of biological macromolecules and argue that RNA 3D structure can be modeled using “protein-like” methodology. We also highlight the areas where the differences between RNA and proteins require the development of RNA-specific solutions.

Figure

Toward a consensus view of duplex RNA flexibility

The structure and flexibility of the RNA duplex has been studied using extended molecular dynamics simulations on four diverse 18-mer oligonucleotides designed to contain many copies of the 10 unique dinucleotide steps in different sequence environments. Simulations were performed using the two most popular force fields for nucleic acids simulations (AMBER and CHARMM) in their latest versions, trying to arrive to a consensus picture of the RNA flexibility. Contrary to what was found for DNA duplex (DNA(2)), no clear convergence is found for the RNA duplex (RNA(2)), but one of the force field seems to agree better with experimental data. MD simulations performed with this force field were used to fully characterize, for the first time to our knowledge, the sequence-dependent elastic properties of RNA duplexes at different levels of resolutions. The flexibility pattern of RNA(2) shows similarities with DNA(2), but also surprising differences, which help us to understand the different biological functions of both molecules. A full mesoscopic model of RNA duplex at different resolution levels is derived to be used for genome-wide description of the flexibility of double-helical fragments of RNA.

Bimolecular Complementation to Visualize Filovirus VP40-Host Complexes in Live Mammalian Cells: Toward the Identification of Budding Inhibitors

Virus-host interactions play key roles in promoting efficient egress of many RNA viruses, including Ebola virus (EBOV or “e”) and Marburg virus (MARV or “m”). Late- (L-) domains conserved in viral matrix proteins recruit specific host proteins, such as Tsg101 and Nedd4, to facilitate the budding process. These interactions serve as attractive targets for the development of broad-spectrum budding inhibitors. A major gap still exists in our understanding of the mechanism of filovirus budding due to the difficulty in detecting virus-host complexes and mapping their trafficking patterns in the natural environment of the cell. To address this gap, we used a bimolecular complementation (BiMC) approach to detect, localize, and follow the trafficking patterns of eVP40-Tsg101 complexes in live mammalian cells. In addition, we used the BiMC approach along with a VLP budding assay to test small molecule inhibitors identified by in silico screening for their ability to block eVP40 PTAP-mediated interactions with Tsg101 and subsequent budding of eVP40 VLPs. We demonstrated the potential broad spectrum activity of a lead candidate inhibitor by demonstrating its ability to block PTAP-dependent binding of HIV-1 Gag to Tsg101 and subsequent egress of HIV-1 Gag VLPs.

Computational methods for studying serpin conformational change and structural plasticity

Currently, over a hundred high-resolution structures of serpins are available, exhibiting a wide range of conformations. However, our understanding of serpin dynamics and conformational change is still limited, mainly due to challenges of monitoring structural changes and characterizing transient conformations using experimental methods. Insight can be provided, however, by employing theoretical and computational approaches. In this chapter, we present an overview of such methods, focusing on molecular dynamics and simulation. As serpin conformational dynamics span a wide range of timescales, we discuss the relative merits of each method and suggest which method is suited to specific conformational phenomena.

Computational studies of colicin insertion into membranes: the closed state

Colicins are water-soluble toxins that, upon interaction with membranes, undergo a conformational change, insert, and form pores in them. Pore formation activity is localized in a bundle of 10 α-helices named the pore-forming domain (PFD). There is evidence that colicins attach to the membrane via a hydrophobic hairpin embedded in the core of the PFD. Two main models have been suggested for the membrane-bound state: penknife and umbrella, differing in regard to the orientation of the hydrophobic hairpin with respect to the membrane. The arrangement of the amphipathic helices has been described as either a compact three-dimensional structure or a two-dimensional array of loosely interacting helices on the membrane surface. Using molecular dynamics simulations with an implicit membrane model, we studied the structure and stability of the conformations proposed earlier for four colicins. We find that colicins are initially driven towards the membrane by electrostatic interactions between basic residues and the negatively charged membrane surface. They do not have a unique binding orientation, but in the predominant orientations the central hydrophobic hairpin is parallel to the membrane. In the inserted state, the estimated free energy tends to be lower for the compact arrangements of the amphipathic helix, but the more expanded ones are in better agreement with experimental distance distributions. The difference in energy between penknife and umbrella conformations is small enough for equilibrium to exist between them. Elongation of the hydrophobic hairpin helices and membrane thinning were found unable to produce stabilization of the transmembrane configuration of the hydrophobic hairpin.

Development of small-molecule PUMA inhibitors for mitigating radiation-induced cell death

PUMA (p53 upregulated modulator of apoptosis) is a Bcl-2 homology 3 (BH3)-only Bcl-2 family member and a key mediator of apoptosis induced by a wide variety of stimuli. PUMA is particularly important in initiating radiation-induced apoptosis and damage in the gastrointestinal and hematopoietic systems. Unlike most BH3-only proteins, PUMA neutralizes all five known antiapoptotic Bcl-2 members though high affinity interactions with its BH3 domain to initiate mitochondria-dependent cell death. Using structural data on the conserved interactions of PUMA with Bcl-2-like proteins, we developed a pharmacophore model that mimics these interactions. In silico screening of the ZINC 8.0 database with this pharmacophore model yielded 142 compounds that could potentially disrupt these interactions. Thirteen structurally diverse compounds with favorable in silico ADME/Toxicity profiles have been retrieved from this set. Extensive testing of these compounds using cell-based and cell-free systems identified lead compounds that confer considerable protection against PUMA-dependent and radiation-induced apoptosis, and inhibit the interaction between PUMA and Bcl-xL.

ROSETTA3: an object-oriented software suite for the simulation and design of macromolecules

We have recently completed a full re-architecturing of the ROSETTA molecular modeling program, generalizing and expanding its existing functionality. The new architecture enables the rapid prototyping of novel protocols by providing easy-to-use interfaces to powerful tools for molecular modeling. The source code of this rearchitecturing has been released as ROSETTA3 and is freely available for academic use. At the time of its release, it contained 470,000 lines of code. Counting currently unpublished protocols at the time of this writing, the source includes 1,285,000 lines. Its rapid growth is a testament to its ease of use. This chapter describes the requirements for our new architecture, justifies the design decisions, sketches out central classes, and highlights a few of the common tasks that the new software can perform.

Modeling of proteins and their assemblies with the integrative modeling platform

To understand the workings of the living cell, we need to characterize protein assemblies that constitute the cell (for example, the ribosome, 26S proteasome, and the nuclear pore complex). A reliable high-resolution structural characterization of these assemblies is frequently beyond the reach of current experimental methods, such as X-ray crystallography, NMR spectroscopy, electron microscopy, footprinting, chemical cross-linking, FRET spectroscopy, small-angle X-ray scattering, and proteomics. However, the information garnered from different methods can be combined and used to build computational models of the assembly structures that are consistent with all of the available datasets. Here, we describe a protocol for this integration, whereby the information is converted to a set of spatial restraints and a variety of optimization procedures can be used to generate models that satisfy the restraints as much as possible. These generated models can then potentially inform about the precision and accuracy of structure determination, the accuracy of the input datasets, and further data generation. We also demonstrate the Integrative Modeling Platform (IMP) software, which provides the necessary computational framework to implement this protocol, and several applications for specific-use cases.

Computational studies of protegrin antimicrobial peptides: a review

Antimicrobial peptides (AMPs) are small, naturally occurring peptides that exhibit strong antibacterial properties generally believed to be a result of selective bacterial membrane disruption. As a result, there has been significant interest in the development of therapeutic antibiotics based on AMPs; however, the poor understanding of the fundamental mechanism of action of these peptides has largely hampered such efforts. We present a summary of computational and theoretical investigations of protegrin, a particularly potent peptide that is both an excellent model for the mechanism of action of AMPs and a promising therapeutic candidate. Experimental investigations have shed light on many of the key steps in the action of protegrin: protegrin monomers are known to dimerize in various lipid environments; protegrin peptides interact strongly with lipid bilayer membranes, particularly anionic lipids; protegrins have been shown to form pores in lipid bilayers, which results in uncontrolled ion transport and may be a key factor in bacterial death. In this work, we present a comprehensive review of the computational and theoretical studies that have complemented and extended the information obtained from experimental work with protegrins, as well as a brief survey of the experimental biophysical studies that are most pertinent to such computational work. We show that a consistent, mechanistic description of the bactericidal mechanism of action of protegrins is emerging, and briefly outline areas where the current understanding is deficient. We hope that the research reviewed herein offers compelling evidence of the benefits of computational investigations of protegrins and other AMPs, as well as providing a useful guide to future work in this area.

Preparation and refinement of model protein-ligand complexes

The formation of ligand-protein complexes are critical for the correct functioning of a cell. The prediction of these interactions is important for our understanding of how the cell works and for the development of new drug molecules. Homology modeling is a method for predicting the structure of a protein based on a crystal structure template. Once a model of the protein is complete, a ligand-docking algorithm predicts the ligand-protein model interaction by searching for the best steric and energetically favorable fit. A refinement of the ligand-binding pocket improves the predicted interactions by considering the flexible nature of the ligand-binding pocket. In this chapter, we describe, from first principles, methods to identify and prepare the ligand-binding pocket in a protein model, to dock the ligand, and refine the resulting complex.

High-performance scalable molecular dynamics simulations of a polarizable force field based on classical Drude oscillators in NAMD

Incorporating the influence of induced polarization in large-scale atomistic molecular dynamics (MD) simulations is a critical challenge in the progress toward computations of increased accuracy. One computationally efficient treatment is based on the classical Drude oscillator in which an auxiliary charged particle is attached by a spring to each nucleus. Here, we report the first implementation of this model in the program NAMD. An extended Lagrangian dynamics with a dual-Langevin thermostat scheme applied to the Drude-nucleus pairs is employed to efficiently generate classical dynamic propagation near the self-consistent field limit. Large-scale MD simulations based on the Drude polarizable force field scale very well on massively distributed supercomputing platforms, the computational demand being only about 50-100% higher than for nonpolarizable models. As an illustration, a large-scale 150 mM NaCl aqueous salt solution is simulated, and the calculated ionic conductivity is shown to be in excellent agreement with experiment.

ROSETTA3: an object-oriented software suite for the simulation and design of macromolecules

We have recently completed a full re-architecturing of the ROSETTA molecular modeling program, generalizing and expanding its existing functionality. The new architecture enables the rapid prototyping of novel protocols by providing easy-to-use interfaces to powerful tools for molecular modeling. The source code of this rearchitecturing has been released as ROSETTA3 and is freely available for academic use. At the time of its release, it contained 470,000 lines of code. Counting currently unpublished protocols at the time of this writing, the source includes 1,285,000 lines. Its rapid growth is a testament to its ease of use. This chapter describes the requirements for our new architecture, justifies the design decisions, sketches out central classes, and highlights a few of the common tasks that the new software can perform.

ROSETTA3: an object-oriented software suite for the simulation and design of macromolecules

We have recently completed a full re-architecturing of the ROSETTA molecular modeling program, generalizing and expanding its existing functionality. The new architecture enables the rapid prototyping of novel protocols by providing easy-to-use interfaces to powerful tools for molecular modeling. The source code of this rearchitecturing has been released as ROSETTA3 and is freely available for academic use. At the time of its release, it contained 470,000 lines of code. Counting currently unpublished protocols at the time of this writing, the source includes 1,285,000 lines. Its rapid growth is a testament to its ease of use. This chapter describes the requirements for our new architecture, justifies the design decisions, sketches out central classes, and highlights a few of the common tasks that the new software can perform.

Trends in template/fragment-free protein structure prediction

Predicting the structure of a protein from its amino acid sequence is a long-standing unsolved problem in computational biology. Its solution would be of both fundamental and practical importance as the gap between the number of known sequences and the number of experimentally solved structures widens rapidly. Currently, the most successful approaches are based on fragment/template reassembly. Lacking progress in template-free structure prediction calls for novel ideas and approaches. This article reviews trends in the development of physical and specific knowledge-based energy functions as well as sampling techniques for fragment-free structure prediction. Recent physical- and knowledge-based studies demonstrated that it is possible to sample and predict highly accurate protein structures without borrowing native fragments from known protein structures. These emerging approaches with fully flexible sampling have the potential to move the field forward.

Computational generation inhibitor-bound conformers of p38 MAP kinase and comparison with experiments

The p38 MAP kinases play a critical role in regulating stress-activated pathways, and serve as molecular targets for controlling inflammatory diseases. Computer-aided efforts for developing p38 inhibitors have been hampered by the necessity to include the enzyme conformational flexibility in ligand docking simulations. A useful strategy in such complicated cases is to perform ensemble-docking provided that a representative set of conformers is available for the target protein either from computations or experiments. We explore here the abilities of two computational approaches, molecular dynamics (MD) simulations and anisotropic network model (ANM) normal mode analysis, for generating potential ligand-bound conformers starting from the apo state of p38, and benchmark them against the space of conformers (or the reference modes of structural changes) inferred from principal component analysis of 134 experimentally resolved p38 kinase structures. ANM-generated conformations are found to provide a significantly better coverage of the inhibitor-bound conformational space observed experimentally, compared to MD simulations performed in explicit water, suggesting that ANM-based sampling of conformations can be advantageously employed as input structural models in docking simulations.

Efficiency of alchemical free energy simulations. II. Improvements for thermodynamic integration

We attempt to optimize the efficiency of thermodynamic integration, as defined by the minimal number of unphysical intermediate states required for the computation of accurate and precise free energy differences. The suitability of various numerical quadrature methods is tested. In particular, we compare the trapezoidal rule, Simpson's rule, Gauss-Legendre, Gauss-Kronrod-Patterson, and Clenshaw-Curtis integration, as well as integration based on a cubic spline approximation of the integrand. We find that Simpson's rule and spline integration are already significantly more efficient that the trapezoidal rule, i.e., correct free energy differences can be obtained using fewer λ-states. We demonstrate that Simpson's rule can be used advantageously with nonequidistant values of the abscissa, which increases the flexibility of the method. Efficiency is enhanced even further if higher order methods, such as Gauss-Legendre, Gauss-Kronrod-Patterson, or Clenshaw-Curtis integration, are used; no more than seven λ-states, which in the case of Clenshaw-Curtis integration include the physical end states, were required for accurate results in all test problems studied. Thus, the performance of thermodynamic integration can equal that of Bennett's acceptance ratio method. We also show, however, that the high efficiency found here relies on the particular functional form of the soft-core potential used; overall, thermodynamic integration is more susceptible to the details of the hybrid Hamiltonian used than Bennett's acceptance ratio method. Therefore, we recommend Bennett's acceptance ratio method as the most robust method to compute alchemical free energy differences; nevertheless, scenarios when thermodynamic integration may be preferable are discussed.

Conformational switching between protein substates studied with 2D IR vibrational echo spectroscopy and molecular dynamics simulations

Myoglobin is an important protein for the study of structure and dynamics. Three conformational substates have been identified for the carbonmonoxy form of myoglobin (MbCO). These are manifested as distinct peaks in the IR absorption spectrum of the CO stretching mode. Ultrafast 2D IR vibrational echo chemical exchange experiments are used to observed switching between two of these substates, A(1) and A(3), on a time scale of <100 ps for two mutants of wild-type Mb. The two mutants are a single mutation of Mb, L29I, and a double mutation, T67R/S92D. Molecular dynamics (MD) simulations are used to model the structural differences between the substates of the two MbCO mutants. The MD simulations are also employed to examine the substate switching in the two mutants as a test of the ability of MD simulations to predict protein dynamics correctly for a system in which there is a well-defined transition over a significant potential barrier between two substates. For one mutant, L29I, the simulations show that translation of the His64 backbone may differentiate the two substates. The simulations accurately reproduce the experimentally observed interconversion time for the L29I mutant. However, MD simulations exploring the same His64 backbone coordinate fail to display substate interconversion for the other mutant, T67R/S92D, thus pointing to the likely complexity of the underlying protein interactions. We anticipate that understanding conformational dynamics in MbCO via ultrafast 2D IR vibrational echo chemical exchange experiments can help to elucidate fast conformational switching processes in other proteins.

Structure of the chromophore binding pocket in the Pr state of plant phytochrome phyA

A homology structural model was generated for plant phytochrome phyA utilizing the crystal structure of the sensory module of cyanobacterial phytochrome Cph1 (Cph1Δ2). As chromophores, either the native phytochromobilin cofactor (PΦB) or phycocyanobilin (PCB), the natural cofactor in Cph1, was incorporated. These homology models were further optimized by molecular dynamics (MD) simulations revealing a satisfying overall agreement with the crystal structure of Cph1Δ2. Notable differences in the PΦB adduct of phyA result from a restructuring of the small helical segment α(7) that leads to displacements of a few amino acids away from the cofactor. This repositioning of residues also include aspartate 218 such that, instead of its carbonyl function as in Cph1Δ2, an additional water molecule forms hydrogen bonds with the ring B and C NH groups. To validate the phyA structural model in the chromophore binding pocket, Raman spectra of the cofactor were calculated by means of the quantum mechanics/molecular mechanics (QM/MM) hybrid methodology and compared with the experimental resonance Raman (RR) spectra. The satisfactory overall agreement between calculated and experimental spectra is taken as an indication for the good quality of the structural model. Moreover, the methine bridge stretching modes and the effects of isotopic labeling at selected positions of the chromophore are very well reproduced to allow confirming even details of the methine bridge geometry as predicted by the homology model. Specifically, it is demonstrated that the experimental RR spectra are consistent with a torsional angle of ring D with respect to ring C that is distinctly higher for phyA-PCB (45°) and phyA-PΦB (42°) than for Cph1Δ2 (30°). Raman spectra calculated from different points of the MD trajectory display variations of the mode frequencies and intensities reflecting the structural fluctuations from snapshot to snapshot. The snapshot spectrum of the lowest energy structure and the sum of all snapshot spectra afford an equally good description of the experimental data. Particularly large variations between the snapshots are noted for the N-H in-plane bending mode of the pyrrole rings B and C, which reflect alterations of the hydrogen bond interactions brought about by fluctuations of water molecules in the cofactor cavity. This overestimation of the water molecule mobility is a consequence of the deficiency of the current QM/MM methodology that, due to the lack of appropriate protein force fields, cannot adequately account for the electrostatics in the cofactor pocket.