CHARMM: the biomolecular simulation program

Tetracycline-tet repressor binding specificity: insights from experiments and simulations

Tetracycline (Tc) antibiotics have been put to new uses in the construction of artificial gene regulation systems, where they bind to the Tet repressor protein (TetR) and modulate its affinity for DNA. Many Tc variants have been produced, both to overcome bacterial resistance and to achieve a broad range of binding strengths. To better understand TetR-Tc binding, we investigate a library of 16 tetracyclines, using fluorescence experiments and molecular dynamics free energy simulations (MDFE). The relative TetR binding free energies are computed by reversibly transforming one Tc variant into another during the simulation, with no adjustable parameters. The chemical variations involve polar and nonpolar substitutions along one entire edge of the elongated Tc structure, which provides many of the protein-ligand contacts. The binding constants span five orders of magnitude. The simulations reproduce the experimental binding free energies, when available, within the uncertainty of either method (+/-0.5 kcal/mol), and reveal many additional details. Contributions of individual Tc substituents are evaluated, along with their additivity and transferability among different positions on the Tc scaffold; differences between D- and B-class repressors are quantified. With increasing computer power, the MDFE approach provides an attractive complement to experiment and should play an increasing role in the understanding and engineering of protein-ligand recognition.

Polarizability rescaling and atom-based Thole scaling in the CHARMM Drude polarizable force field for ethers

Within the CHARMM polarizable force field based on the classical Drude oscillator, atomic polarizabilities are derived via fitting to ab initio calculated data on isolated gas phase molecules, with an empirical scaling factor applied to account for differences between the gas and condensed phases. In the development of polarizable models for the ethers, a polarizability scaling factor of 0.7 was previously applied [Vorobyov et al. J Comput Chem 3:1120-1133, 2007]. While the resulting force field models gave good agreement with a variety of experimental data, they systematically underestimated the liquid phase dielectric constants. Here, a new CHARMM polarizable model is developed for the ethers, employing a polarizability scaling factor of 0.85 and including atom-based Thole scale factors recently introduced into the CHARMM Drude polarizable force field [Harder et al. J Phys Chem B 112:3509-3521, 2008]. The new model offers a significant improvement in the reproduction of liquid phase dielectric constants, while maintaining the good agreement of the previous model with all other experimental and quantum mechanical data, highlighting the sensitivity of liquid phase properties to the choice of atomic polarizability parameters.

Accurate Calculation of Hydration Free Energies using Pair-Specific Lennard-Jones Parameters in the CHARMM Drude Polarizable Force Field

Lennard-Jones (LJ) parameters for a variety of model compounds have previously been optimized within the CHARMM Drude polarizable force field to reproduce accurately pure liquid phase thermodynamic properties as well as additional target data. While the polarizable force field resulting from this optimization procedure has been shown to satisfactorily reproduce a wide range of experimental reference data across numerous series of small molecules, a slight but systematic overestimate of the hydration free energies has also been noted. Here, the reproduction of experimental hydration free energies is greatly improved by the introduction of pair-specific LJ parameters between solute heavy atoms and water oxygen atoms that override the standard LJ parameters obtained from combining rules. The changes are small and a systematic protocol is developed for the optimization of pair-specific LJ parameters and applied to the development of pair-specific LJ parameters for alkanes, alcohols and ethers. The resulting parameters not only yield hydration free energies in good agreement with experimental values, but also provide a framework upon which other pair-specific LJ parameters can be added as new compounds are parametrized within the CHARMM Drude polarizable force field. Detailed analysis of the contributions to the hydration free energies reveals that the dispersion interaction is the main source of the systematic errors in the hydration free energies. This information suggests that the systematic error may result from problems with the LJ combining rules and is combined with analysis of the pair-specific LJ parameters obtained in this work to identify a preliminary improved combining rule.

RNA polymerase II flexibility during translocation from normal mode analysis

The structural dynamics in eukaryotic RNA polymerase II (RNAPII) is described from computational normal mode analysis based on a series of crystal structures of pre- and post-translocated states with open and closed trigger loops. Conserved modes are identified that involve translocation of the nucleic acid complex coupled to motions of the enzyme, in particular in the clamp and jaw domains of RNAPII. A combination of these modes is hypothesized to be involved during active transcription. The NMA modes indicate furthermore that downstream DNA translocation may occur separately from DNA:RNA hybrid translocation. A comparison of the modes between different states of RNAPII suggests that productive translocation requires an open trigger loop and is inhibited by the presence of an NTP in the active site. This conclusion is also supported by a comparison of the overall flexibility in terms of root mean square fluctuations.

Preferential interactions of guanidinum ions with aromatic groups over aliphatic groups

Small angle neutron scattering (SANS) and molecular dynamics (MD) simulations were used to characterize the long-range structuring (aggregation) of aqueous solutions of isopropanol (IPA) and pyridine and the effect on structuring of guanidinium chloride (GdmCl). These solutes serve as highly soluble analogs of the nonpolar aliphatic (IPA) and aromatic (pyridine) side chains of proteins. SANS data showed that isopropanol and pyridine both form clusters in water resulting from interaction between nonpolar groups of the solutes, with pyridine aggregation producing longer-range structuring than isopropanol in 3 m solutions. Addition of GdmCl at 3 m concentration considerably reduced pyridine aggregation but had no effect on isopropanol aggregation. MD simulations of these solutions support the conclusion that long-range structuring involves hydrophobic solute interactions and that Gdm(+) interacts with the planar pyridine group to suppress pyridine-pyridine interactions in solution. Hydrophobic interactions involving the aliphatic groups of isopropanol were unaffected by GdmCl, indicating that the planar and weakly hydrated Gdm(+) cation cannot make productive interactions with the highly curved or "lumpy" aliphatic groups of this solute. These observations support the conclusion that the effects of Gdm(+) ions on protein-stabilizing interactions involving aromatic amino acid side chains make significant contributions to the denaturant activity of GdmCl, whereas interactions with the "lumpy" aliphatic side chains are likely to be less important.

A conserved protonation-induced switch can trigger "ionic-lock" formation in adrenergic receptors

The mechanism of signal transduction in G-protein-coupled receptors (GPCRs) is a crucial step in cell signaling. However, the molecular details of this process are still largely undetermined. Carrying out submicrosecond molecular dynamics simulations of beta-adrenergic receptors, we found that cooperation between a number of highly conserved residues is crucial to alter the equilibrium between the active state and the inactive state of diffusible ligand GPCRs. In particular, "ionic-lock" formation in beta-adrenergic receptors is directly correlated with the protonation state of a highly conserved aspartic acid residue [Asp(2.50)] even though the two sites are located more than 20 A away from each other. Internal polar residues, acting as local microswitches, cooperate to propagate the signal from Asp(2.50) to the G-protein interaction site at the helix III-helix VI interface. Evolutionarily conserved differences between opsin and non-opsin GPCRs in the surrounding of Asp(2.50) influence the acidity of this residue and can thus help in rationalizing the differences in constitutive activity of class A GPCRs.

The relationship between curvature, flexibility and persistence length in the tropomyosin coiled-coil

The inherent flexibility of rod-like tropomyosin coiled-coils is a significant factor that constrains tropomyosin's complex positional dynamics on actin filaments. Flexibility of elongated straight molecules typically is assessed by persistence length, a measure of lengthwise thermal bending fluctuations. However, if a molecule's equilibrium conformation is curved, this formulation yields an "apparent" persistence length ( approximately 100nm for tropomyosin), measuring deviations from idealized straight conformations which then overestimate actual dynamic flexibility. To obtain the "dynamic" persistence length, a true measurement of flexural stiffness, the average curvature of the molecule must be taken into account. Different methods used in our studies for measuring the dynamic persistence length directly from Molecular Dynamics (MD) simulations of tropomyosin are described here in detail. The dynamic persistence length found, 460+/-40nm, is approximately 12-times longer than tropomyosin and 5-times the apparent persistence length, showing that tropomyosin is considerably stiffer than previously thought. The longitudinal twisting behavior of tropomyosin during MD shows that the amplitude of end-to-end twisting fluctuation is approximately 30 degrees when tropomyosin adopts its near-average conformation. The measured bending and twisting flexibilities are used to evaluate different models of tropomyosin motion on F-actin.

Evaporation of solvent molecules by ultrafast heating: effect on conformation of solvated protein

Using molecular dynamics simulation, we compare two cases of ultrafast heating of a small water droplet containing a solvated protein (echistatin). If the water temperature after irradiation is above the critical temperature, explosive boiling liberates the protein within some 10 ps of its hydration shell, while its temperature remains relatively low. By comparing with the case where the water shell is heated to the same final temperature, but without complete evaporation, we demonstrate that the protein conformation is governed by the hydration shell rather than by the protein temperature.

Internal proton transfer in the external pyridoxal 5'-phosphate Schiff base in dopa decarboxylase

Combined quantum mechanical and molecular mechanical (QM/MM) simulations of dopa decarboxylase have been carried out to elucidate the factors that contribute to the tautomeric equilibrium of the intramolecular proton transfer in the external PLP-L-dopa Schiff base. The presence of a carboxylate anion on the alpha-carbon of the Schiff base stabilizes the zwitterions and shifts the equilibrium in favor of the oxoenamine tautomer (protonated Schiff base). Moreover, protonation of the PLP pyridine nitrogen further drives the equilibrium toward the oxoenamine direction. On the other hand, solvent effects favor the hydroxyimine configuration, although the equilibrium favors the oxoenamine isomer with a methyl group as the substituent on the imino nitrogen. In dopa decarboxylase, the hydroxyimine form of the PLP(H+)-L-dopa Schiff base is predicted to be the major isomer with a relative free energy of -1.3 kcal/mol over that of the oxoenamine isomer. Both Asp271 and Lys303 stabilize the hydroxyimine configuration through hydrogen-bonding interactions with the pyridine nitrogen of the PLP and the imino nitrogen of the Schiff base, respectively. Interestingly, Thr246 plays a double role in the intramolecular proton transfer process, in which it initially donates a hydrogen bond to the phenolate oxygen in the oxoenamine configuration and then switches to a hydrogen bond acceptor from the phenolic hydroxyl group in the hydroxyimine tautomer.

Crystal structure of the halotolerant gamma-glutamyltranspeptidase from Bacillus subtilis in complex with glutamate reveals a unique architecture of the solvent-exposed catalytic pocket

gamma-Glutamyltranspeptidase (GGT; EC 2.3.2.2), an enzyme found in organisms from bacteria to mammals and plants, plays a central role in glutathione metabolism. Structural studies of GGTs from Escherichia coli and Helicobacter pylori have revealed detailed molecular mechanisms of catalysis and maturation. In these two GGTs, highly conserved residues form the catalytic pockets, conferring the ability of the loop segment to shield the bound gamma-glutamyl moiety from the solvent. Here, we have examined the Bacillus subtilis GGT, which apparently lacks the amino acids corresponding to the lid-loop that are present in mammalian and plant GGTs as well as in most bacterial GGTs. Another remarkable feature of B. subtilis GGT is its salt tolerance; it retains 86% of its activity even in 3 m NaCl. To better understand these characteristics of B. subtilis GGT, we determined its crystal structure in complex with glutamate, a product of the enzymatic reaction, at 1.95 A resolution. This structure revealed that, unlike the E. coli and H. pylori GGTs, the catalytic pocket of B. subtilis GGT has no segment that covers the bound glutamate; consequently, the glutamate is exposed to solvent. Furthermore, calculation of the electrostatic potential showed that strong acidic patches were distributed on the surface of the B. subtilis GGT, even under high-salt conditions, and this may allow the protein to remain in the hydrated state and avoid self-aggregation. The structural findings presented here have implications for the molecular mechanism of GGT.

Molecular dynamics simulations of nonpolarizable inorganic salt solution interfaces: NaCl, NaBr, and NaI in transferable intermolecular potential 4-point with charge dependent polarizability (TIP4P-QDP) water

We present molecular dynamics simulations of the liquid-vapor interface of 1M salt solutions of nonpolarizable NaCl, NaBr, and NaI in polarizable transferable intermolecular potential 4-point with charge dependent polarizability water [B. A. Bauer et al., J. Chem. Theory Comput. 5, 359 (2009)]; this water model accommodates increased solvent polarizability (relative to the condensed phase) in the interfacial and vapor regions. We employ fixed-charge ion models developed in conjunction with the TIP4P-QDP water model to reproduce ab initio ion-water binding energies and ion-water distances for isolated ion-water pairs. The transferability of these ion models to the condensed phase was validated with hydration free energies computed using thermodynamic integration (TI) and appropriate energy corrections. Density profiles of Cl(-), Br(-), and I(-) exhibit charge layering in the interfacial region; anions and cation interfacial probabilities show marked localization, with the anions penetrating further toward the vapor than the cations. Importantly, in none of the cases studied do anions favor the outermost regions of the interface; there is always an aqueous region between the anions and vapor phase. Observed interfacial charge layering is independent of the strength of anion-cation interactions as manifest in anion-cation contact ion pair peaks and solvent separated ion pair peaks; by artificially modulating the strength of anion-cation interactions (independent of their interactions with solvent), we find little dependence on charge layering particularly for the larger iodide anion. The present results reiterate the widely held view of the importance of solvent and ion polarizability in mediating specific anion surface segregation effects. Moreover, due to the higher parametrized polarizability of the TIP4P-QDP condensed phase {1.31 A(3) for TIP4P-QDP versus 1.1 A(3) (TIP4P-FQ) and 0.87 A(3) (POL3) [Ponder and Case, Adv. Protein Chem. 66, 27 (2003)]} based on ab initio calculations of the condensed-phase polarizability reduction in liquid water, the present simulations highlight the role of water polarizability in inducing water molecular dipole moments parallel to the interface normal (and within the interfacial region) so as to favorably oppose the macrodipole generated by the separation of anion and cation charge. Since the TIP4P-QDP water polarizability approaches that of the experimental vapor phase value for water, the present results suggest a fundamental role of solvent polarizability in accommodating the large spatial dipole generated by the separation of ion charges. The present results draw further attention to the question of what exact value of condensed phase water polarizability to incorporate in classical polarizable water force fields.

Structure and hydration of membranes embedded with voltage-sensing domains

Despite the growing number of atomic-resolution membrane protein structures, direct structural information about proteins in their native membrane environment is scarce. This problem is particularly relevant in the case of the highly charged S1-S4 voltage-sensing domains responsible for nerve impulses, where interactions with the lipid bilayer are critical for the function of voltage-activated ion channels. Here we use neutron diffraction, solid-state nuclear magnetic resonance (NMR) spectroscopy and molecular dynamics simulations to investigate the structure and hydration of bilayer membranes containing S1-S4 voltage-sensing domains. Our results show that voltage sensors adopt transmembrane orientations and cause a modest reshaping of the surrounding lipid bilayer, and that water molecules intimately interact with the protein within the membrane. These structural findings indicate that voltage sensors have evolved to interact with the lipid membrane while keeping energetic and structural perturbations to a minimum, and that water penetrates the membrane, to hydrate charged residues and shape the transmembrane electric field.

Global dynamics of proteins: bridging between structure and function

Biomolecular systems possess unique, structure-encoded dynamic properties that underlie their biological functions. Recent studies indicate that these dynamic properties are determined to a large extent by the topology of native contacts. In recent years, elastic network models used in conjunction with normal mode analyses have proven to be useful for elucidating the collective dynamics intrinsically accessible under native state conditions, including in particular the global modes of motions that are robustly defined by the overall architecture. With increasing availability of structural data for well-studied proteins in different forms (liganded, complexed, or free), there is increasing evidence in support of the correspondence between functional changes in structures observed in experiments and the global motions predicted by these coarse-grained analyses. These observed correlations suggest that computational methods may be advantageously employed for assessing functional changes in structure and allosteric mechanisms intrinsically favored by the native fold.

Coordinate-dependent diffusion in protein folding

Diffusion on a low-dimensional free-energy surface is a remarkably successful model for the folding dynamics of small single-domain proteins. Complicating the interpretation of both simulations and experiments is the expectation that the effective diffusion coefficient D will in general depend on the position along the folding coordinate, and this dependence may vary for different coordinates. Here we explore the position dependence of D, its connection to protein internal friction, and the consequences for the interpretation of single-molecule experiments. We find a large decrease in D from unfolded to folded, for reaction coordinates that directly measure fluctuations in Cartesian configuration space, including those probed in single-molecule experiments. In contrast, D is almost independent of Q, the fraction of native amino acid contacts: Near the folded state, small fluctuations in position cause large fluctuations in Q, and vice versa for the unfolded state. In general, position-dependent free energies and diffusion coefficients for any two good reaction coordinates that separate reactant, product, and transition states, are related by a simple transformation, as we demonstrate. With this transformation, we obtain reaction coordinates with position-invariant D. The corresponding free-energy surfaces allow us to justify the assumptions used in estimating the speed limit for protein folding from experimental measurements of the reconfiguration time in the unfolded state, and also reveal intermediates hidden in the original free-energy projection. Lastly, we comment on the design of future single-molecule experiments that probe the position dependence of D directly.

Role of the adenine ligand on the stabilization of the secondary and tertiary interactions in the adenine riboswitch

Riboswitches are RNA-based genetic control elements that function via a conformational transition mechanism when a specific target molecule binds to its binding pocket. To facilitate an atomic detail interpretation of experimental investigations on the role of the adenine ligand on the conformational properties and kinetics of folding of the add adenine riboswitch, we performed molecular dynamics simulations in both the presence and the absence of the ligand. In the absence of ligand, structural deviations were observed in the J23 junction and the P1 stem. Destabilization of the P1 stem in the absence of ligand involves the loss of direct stabilizing interactions with the ligand, with additional contributions from the J23 junction region. The J23 junction of the riboswitch is found to be more flexible, and the tertiary contacts among the junction regions are altered in the absence of the adenine ligand; results suggest that the adenine ligand associates and dissociates from the riboswitch in the vicinity of J23. Good agreement was obtained with the experimental data with the results indicating dynamic behavior of the adenine ligand on the nanosecond time scale to be associated with the dynamic behavior of hydrogen bonding with the riboswitch. Results also predict that direct interactions of the adenine ligand with U74 of the riboswitch are not essential for stable binding although it is crucial for its recognition. The possibility of methodological artifacts and force-field inaccuracies impacting the present observations was checked by additional molecular dynamics simulations in the presence of 2,6-diaminopurine and in the crystal environment.

Kinetic characterization of xenobiotic reductase A from Pseudomonas putida 86

Xenobiotic reductase A (XenA) from Pseudomonas putida is a member of the old-yellow-enzyme family of flavin-containing enzymes and catalyzes the NADH/NADPH-dependent reduction of various substrates, including 8-hydroxycoumarin and 2-cyclohexenone. Here we present a kinetic and thermodynamic analysis of XenA. In the reductive half-reaction, complexes of oxidized XenA with NADH or NADPH form charge-transfer (CT) intermediates with increased absorption around 520-560 nm, which occurs with a second-order rate constant of 9.4 x 10(5) M(-1) s(-1) with NADH and 6.4 x 10(5) M(-1) s(-1) with NADPH, while its disappearance is controlled by a rate constant of 210-250 s(-1) with both substrates. Transfer of hydride from NADPH proceeds 24 times more rapidly than from NADH. This modest kinetic preference of XenA for NADPH is unlike the typical discrimination between NADH and NADPH by binding affinity. Docking studies combined with electrostatic energy calculations indicate that the 2'-phosphate group attached to the adenine moiety of NADPH is responsible for this difference. The reductions of 2-cyclohexenone and coumarin in the oxidative half-reaction are both concentration-dependent under the assay conditions and reveal a more than 50-fold larger limiting rate constant for the reduction of 2-cyclohexenone compared to that of coumarin. Our work corroborates the link between XenA and other members of the old-yellow-enzyme family but demonstrates several differences in the reactivity of these enzymes.

Protein-protein docking using region-based 3D Zernike descriptors

Background

Protein-protein interactions are a pivotal component of many biological processes and mediate a variety of functions. Knowing the tertiary structure of a protein complex is therefore essential for understanding the interaction mechanism. However, experimental techniques to solve the structure of the complex are often found to be difficult. To this end, computational protein-protein docking approaches can provide a useful alternative to address this issue. Prediction of docking conformations relies on methods that effectively capture shape features of the participating proteins while giving due consideration to conformational changes that may occur.

Results

We present a novel protein docking algorithm based on the use of 3D Zernike descriptors as regional features of molecular shape. The key motivation of using these descriptors is their invariance to transformation, in addition to a compact representation of local surface shape characteristics. Docking decoys are generated using geometric hashing, which are then ranked by a scoring function that incorporates a buried surface area and a novel geometric complementarity term based on normals associated with the 3D Zernike shape description. Our docking algorithm was tested on both bound and unbound cases in the ZDOCK benchmark 2.0 dataset. In 74% of the bound docking predictions, our method was able to find a near-native solution (interface C-αRMSD ≤ 2.5 Å) within the top 1000 ranks. For unbound docking, among the 60 complexes for which our algorithm returned at least one hit, 60% of the cases were ranked within the top 2000. Comparison with existing shape-based docking algorithms shows that our method has a better performance than the others in unbound docking while remaining competitive for bound docking cases.

Conclusion

We show for the first time that the 3D Zernike descriptors are adept in capturing shape complementarity at the protein-protein interface and useful for protein docking prediction. Rigorous benchmark studies show that our docking approach has a superior performance compared to existing methods.

Structure-based optimization of potent and selective inhibitors of the tyrosine kinase erythropoietin producing human hepatocellular carcinoma receptor B4 (EphB4)

The tyrosine kinase EphB4 is an attractive target for drug design because of its recognized role in cancer-related angiogenesis. Recently, a series of commercially available xanthine derivatives were identified as micromolar inhibitors of EphB4 by high-throughput fragment-based docking into the ATP-binding site of the kinase domain. Here, we have exploited the binding mode obtained by automatic docking for the optimization of these EphB4 inhibitors by chemical synthesis. Addition of only two heavy atoms, methyl and hydroxyl groups, to compound 4 has yielded the single-digit nanomolar inhibitor 66, with a remarkable improvement of the ligand efficiency from 0.26 to 0.37 kcal/(mol per non-hydrogen atom). Compound 66 shows very high affinity for a few other tyrosine kinases with threonine as gatekeeper residue (Abl, Lck, and Src). On the other hand, it is selective against kinases with a larger gatekeeper. A 45 ns molecular dynamics (MD) simulation of the complex of EphB4 and compound 66 provides further validation of the binding mode obtained by fragment-based docking.

Differential quantum tunneling contributions in nitroalkane oxidase catalyzed and the uncatalyzed proton transfer reaction

The proton transfer reaction between the substrate nitroethane and Asp-402 catalyzed by nitroalkane oxidase and the uncatalyzed process in water have been investigated using a path-integral free-energy perturbation method. Although the dominating effect in rate acceleration by the enzyme is the lowering of the quasiclassical free energy barrier, nuclear quantum effects also contribute to catalysis in nitroalkane oxidase. In particular, the overall nuclear quantum effects have greater contributions to lowering the classical barrier in the enzyme, and there is a larger difference in quantum effects between proton and deuteron transfer for the enzymatic reaction than that in water. Both experiment and computation show that primary KIEs are enhanced in the enzyme, and the computed Swain-Schaad exponent for the enzymatic reaction is exacerbated relative to that in the absence of the enzyme. In addition, the computed tunneling transmission coefficient is approximately three times greater for the enzyme reaction than the uncatalyzed reaction, and the origin of the difference may be attributed to a narrowing effect in the effective potentials for tunneling in the enzyme than that in aqueous solution.

Entrapment and structure of an extrahelical guanine attempting to enter the active site of a bacterial DNA glycosylase, MutM

MutM, a bacterial DNA glycosylase, protects genome integrity by catalyzing glycosidic bond cleavage of 8-oxoguanine (oxoG) lesions, thereby initiating base excision DNA repair. The process of searching for and locating oxoG lesions is especially challenging, because of the close structural resemblance of oxoG to its million-fold more abundant progenitor, G. Extrusion of the target nucleobase from the DNA double helix to an extrahelical position is an essential step in lesion recognition and catalysis by MutM. Although the interactions between the extruded oxoG and the active site of MutM have been well characterized, little is known in structural detail regarding the interrogation of extruded normal DNA bases by MutM. Here we report the capture and structural elucidation of a complex in which MutM is attempting to present an undamaged G to its active site. The structure of this MutM-extrahelical G complex provides insights into the mechanism MutM employs to discriminate against extrahelical normal DNA bases and into the base extrusion process in general.

N-acylpolyamine inhibitors of HDM2 and HDMX binding to p53

Selective inhibition of protein-protein interactions important for cellular processes could lead to the development of new therapies against disease. In the area of cancer, overexpression of the proteins human double minute 2 (HDM2) and its homolog HDMX has been linked to tumor aggressiveness. Both HDM2 and HDMX bind to p53 and prevent cell cycle arrest or apoptosis in damaged cells. Developing a strategy to simultaneously prevent the binding of both HDM2 and HDMX to p53 is an essential feature of inhibitors to restore p53 activity in a number of different cancers. Inhibition of protein-protein interactions with synthetic molecules is an emerging area of research that requires new inhibitors tailored to mimic the types of interfaces between proteins. Our strategy to create inhibitors of protein-protein interactions is to develop a non-natural scaffold that may be used as a starting point to identify important molecular components necessary for inhibition. In this study, we report an N-acylpolyamine (NAPA) scaffold that supports numerous sidechains in a compact atomic arrangement. NAPAs were constructed by a series of reductive aminations between amino acid derivatives followed by acylation at the resulting secondary amine. An optimized NAPA was able to equally inhibit the association of both HDM2 and HDMX with p53. Our results demonstrate some of the challenges associated with targeting multiple protein-protein interactions involved in overlapping cellular processes.

Predicted structural basis for CD1c presentation of mycobacterial branched polyketides and long lipopeptide antigens

CD1 proteins mediate the trafficking and presentation of a diverse range of lipid antigens to T-cell receptors, and thus play a key role in our adaptive immune system. Crystal structures of several CD1 isoforms reveal a highly conserved tertiary structure, but also great variability in the anatomy of their binding pockets, reflecting their distinct ligand specificity. The structure of one important member of the family, CD1c, remains unknown. CD1c is of great interest as it can present an unusual and potent lipid antigen, mannosyl-beta(1)-phosphomycoketide (MPM) from Mycobacterium tuberculosis, the causative agent of tuberculosis. CD1c has also been reported to present acetylated 12-amino-acid-long peptides (lipo-12), an observation with broad immunological implications but difficult to rationalize on structural grounds. To gain insights into the structural basis for the ligand specificity of CD1c, we have generated an atomic model of its binding domain using a detailed position-specific multiple-template homology modeling approach. This model reveals structural features unique to this isoform, particularly with regard to the so-called pocket F', which provide a compelling rationale for the ability of CD1c to bind not only branched alkyl chains such as in MPM, but also long lipopeptides comparable to those presented by MHC proteins. A model of CD1c with bound MPM was constructed and analyzed through molecular dynamics simulations, showing marked structural stability in the time-scale of 100 ns. A model of CD1c in complex with lipo-12 is also presented.

Flaviviral protease inhibitors identified by fragment-based library docking into a structure generated by molecular dynamics

Fragment-based docking was used to select a conformation for virtual screening from a molecular dynamics trajectory of the West Nile virus nonstructural 3 protease. This conformation was chosen from an ensemble of 100 molecular dynamics snapshots because it optimally accommodates benzene, the most common ring in known drugs, and two positively charged fragments (methylguanidinium and 2-phenylimidazoline). The latter fragments were used as probes because of the large number of hydrogen bond acceptors in the substrate binding site of the protease. Upon high-throughput docking of a diversity set of 18,694 molecules and pose filtering, only five compounds were chosen for experimental validation, and two of them are active in the low micromolar range in an enzymatic assay and a tryptophan fluorescence quenching assay. Evidence for specific binding to the protease active site is provided by nuclear magnetic resonance spectroscopy. The two inhibitors have different scaffolds (diphenylurea and diphenyl ester) and are promising lead candidates because they have a molecular weight of about 300 Da.

CHARMM additive all-atom force field for aldopentofuranoses, methyl-aldopentofuranosides, and fructofuranose

An additive all-atom empirical force field for aldopentofuranoses, methyl-aldopentofuranosides (Me-aldopentofuranosides), and fructofuranose carbohydrates, compatible with existing CHARMM carbohydrate parameters, is presented. Building on existing parameters transferred from cyclic ethers and hexopyranoses, parameters were further developed using target data for complete furanose carbohydrates as well as O-methyl tetrahydrofuran. The bond and angle equilibrium parameters were adjusted to reproduce target geometries from a survey of furanose crystal structures, and dihedral parameters were fit to over 1700 quantum mechanical (QM) MP2/cc-pVTZ//MP2/6-31G(d) conformational energies. The conformational energies were for a variety of complete furanose monosaccharides and included two-dimensional ring pucker energy surfaces. Bonded parameter optimization led to the correct description of the ring pucker for a large set of furanose compounds, while furanose-water interaction energies and distances reproduced QM HF/6-31G(d) results for a number of furanose monosaccharides, thereby validating the nonbonded parameters. Crystal lattice unit cell parameters and volumes, aqueous-phase densities, and aqueous NMR ring pucker and exocyclic data were used to validate the parameters in condensed-phase environments. Conformational sampling analysis of the ring pucker and exocyclic group showed excellent agreement with experimental NMR data, demonstrating that the conformational energetics in aqueous solution are accurately described by the optimized force field. Overall, the parameters reproduce available experimental data well and are anticipated to be of utility in future computational studies of carbohydrates, including in the context of proteins, nucleic acids, and/or lipids when combined with existing CHARMM biomolecular force fields.

How does a simplified-sequence protein fold?

To investigate a putatively primordial protein we have simplified the sequence of a 56-residue alpha/beta fold (the immunoglobulin-binding domain of protein G) by replacing it with polyalanine, polythreonine, and diglycine segments at regions of the sequence that in the folded structure are alpha-helical, beta-strand, and turns, respectively. Remarkably, multiple folding and unfolding events are observed in a 15-micros molecular dynamics simulation at 330 K. The most stable state (populated at approximately 20%) of the simplified-sequence variant of protein G has the same alpha/beta topology as the wild-type but shows the characteristics of a molten globule, i.e., loose contacts among side chains and lack of a specific hydrophobic core. The unfolded state is heterogeneous and includes a variety of alpha/beta topologies but also fully alpha-helical and fully beta-sheet structures. Transitions within the denatured state are very fast, and the molten-globule state is reached in <1 micros by a framework mechanism of folding with multiple pathways. The native structure of the wild-type is more rigid than the molten-globule conformation of the simplified-sequence variant. The difference in structural stability and the very fast folding of the simplified protein suggest that evolution has enriched the primordial alphabet of amino acids mainly to optimize protein function by stabilization of a unique structure with specific tertiary interactions.

Free energy and kinetics of conformational transitions from Voronoi tessellated milestoning with restraining potentials

Milestoning is a method aimed at reconstructing the statistical properties of the long-time dynamics of a system by exploiting the crossing statistics of a set of hypersurfaces, called the "milestones", placed along the reaction coordinate [Faradjian and Elber, J. Chem. Phys.2004, 120, 10880]. Recently, Vanden-Eijnden and Venturoli [J. Chem. Phys.2009, 130, 194101] showed that when a complete Voronoi tessellation of the configurational space is available, milestoning can be reformulated to utilise the statistics from a series of independent simulations, each confined within a single cell via strict reflections at the boundaries. As a byproduct, this "Voronoi tessellated milestoning" method also permits to compute the free energy of the tessellation. Here, the method is extended to support the usage of differentiable restraining potentials to confine the trajectories within each cell.

CHARMM Additive All-Atom Force Field for Glycosidic Linkages between Hexopyranoses

We present an extension of the CHARMM hexopyranose monosaccharide additive all-atom force field to enable modeling of glycosidic-linked hexopyranose polysaccharides. The new force field parameters encompass 1→1, 1→2, 1→3, 1→4, and 1→6 hexopyranose glycosidic linkages, as well as O-methylation at the C(1) anomeric carbon, and are developed to be consistent with the CHARMM all-atom biomolecular force fields for proteins, nucleic acids, and lipids. The parameters are developed in a hierarchical fashion using model compounds containing the key atoms in the full carbohydrates, in particular O-methyl-tetrahydropyran and glycosidic-linked dimers consisting of two molecules of tetrahyropyran or one of tetrahydropyran and one of cyclohexane. Target data for parameter optimization include full two-dimensional energy surfaces defined by the Φ/Ψ glycosidic dihedral angles in the disaccharide analogs as determined by quantum mechanical MP2/cc-pVTZ single point energies on MP2/6-31G(d) optimized structures (MP2/cc-pVTZ//MP2/6-31G(d)). In order to achieve balanced, transferable dihedral parameters for the Φ/Ψ glycosidic dihedral angles, surfaces for all possible chiralities at the ring carbon atoms involved in the glycosidic linkages are considered, resulting in over 5000 MP2/cc-pVTZ//MP2/6-31G(d) conformational energies. Also included as target data are vibrational frequencies, pair interaction energies and distances with water molecules, and intramolecular geometries including distortion of the glycosidic valence angle as a function of the glycosidic dihedral angles. The model-compound optimized force field parameters are validated on full disaccharides through comparison of molecular dynamics results to available experimental data. Good agreement is achieved with experiment for a variety of properties including crystal cell parameters and intramolecular geometries, aqueous densities, and aqueous NMR coupling constants associated with the glycosidic linkage. The newly-developed parameters allow for the modeling of linear, branched, and cyclic hexopyranose glycosides both alone and in heterogenous systems including proteins, nucleic acids and/or lipids when combined with existing CHARMM biomolecular force fields.

Mechanical signaling on the single protein level studied using steered molecular dynamics

Efficient communication between the cell and its external environment is of the utmost importance to the function of multicellular organisms. While signaling events can be generally characterized as information exchange by means of controlled energy conversion, research efforts have hitherto mainly been concerned with mechanisms involving chemical and electrical energy transfer. Here, we review recent computational efforts addressing the function of mechanical force in signal transduction. Specifically, we focus on the role of steered molecular dynamics (SMD) simulations in providing details at the atomic level on a group of protein domains, which play a fundamental role in signal exchange by responding properly to mechanical strain. We start by giving a brief introduction to the SMD technique and general properties of mechanically stable protein folds, followed by specific examples illustrating three general regimes of signal transfer utilizing mechanical energy: purely mechanical, mechanical to chemical, and chemical to mechanical. Whenever possible the physiological importance of the example at hand is stressed to highlight the diversity of the processes in which mechanical signaling plays a key role. We also provide an overview of future challenges and perspectives for this rapidly developing field.

Complex formation between calmodulin and a peptide from the intracellular loop of the gap junction protein connexin43: Molecular conformation and energetics of binding

Gap junctions are formed by a family of transmembrane proteins, connexins. Connexin43 is a widely studied member of the family, being ubiquitously expressed in a variety of tissues and a target of a large number of disease mutations. The intracellular loop of connexin43 has been shown to include a calmodulin binding domain, but detailed 3-dimensional data on the structure of the complex are not available. In this study, we used a synthetic peptide from this domain to reveal the conformation of the calmodulin-peptide complex by small angle X-ray scattering. Upon peptide binding, calmodulin lost its dumbbell shape, adopting a more globular conformation. We also studied the energetics of the interaction using calorimetry and computational methods. All our data indicate that calmodulin binds to the peptide from cx43 in the classical 'collapsed' conformation.

Isotropic periodic sum of electrostatic interactions for polar systems

Isotropic periodic sum (IPS) is a method to calculate long-range interactions based on homogeneity of simulation systems. Long-range interactions are represented by interactions with isotropic periodic images of a defined local region and can be reduced to short ranged IPS potentials. The original electrostatic three-dimensional (3D)-IPS potential was derived based on a nonpolar homogeneous approximation and its application is limited to nonpolar or weak polar systems. This work derived a polar electrostatic 3D-IPS potential based on polar interactions. For the convenience of application, polynomial functions with rationalized coefficients are proposed for electrostatic and Lennard-Jones 3D-IPS potentials. Model systems of various polarities and several commonly used solvent systems are simulated to evaluate the 3D-IPS potentials. It is demonstrated that for polar systems the polar electrostatic 3D-IPS potential has much improved accuracy as compared to the nonpolar 3D-IPS potential. For homogeneous systems, the polar electrostatic 3D-IPS potential with a local region radius or cutoff distance of as small as 10 A can satisfactorily reproduce energetic, structural, and dynamic properties from the particle-meshed-Ewald method. For both homogeneous and heterogeneous systems, the 3D-IPS/discrete fast Fourier transform method using either the nonpolar or the polar electrostatic 3D-IPS potentials results in very similar simulation results.

Analysis of the free-energy surface of proteins from reversible folding simulations

Computer generated trajectories can, in principle, reveal the folding pathways of a protein at atomic resolution and possibly suggest general and simple rules for predicting the folded structure of a given sequence. While such reversible folding trajectories can only be determined ab initio using all-atom transferable force-fields for a few small proteins, they can be determined for a large number of proteins using coarse-grained and structure-based force-fields, in which a known folded structure is by construction the absolute energy and free-energy minimum. Here we use a model of the fast folding helical λ-repressor protein to generate trajectories in which native and non-native states are in equilibrium and transitions are accurately sampled. Yet, representation of the free-energy surface, which underlies the thermodynamic and dynamic properties of the protein model, from such a trajectory remains a challenge. Projections over one or a small number of arbitrarily chosen progress variables often hide the most important features of such surfaces. The results unequivocally show that an unprojected representation of the free-energy surface provides important and unbiased information and allows a simple and meaningful description of many-dimensional, heterogeneous trajectories, providing new insight into the possible mechanisms of fast-folding proteins.

The process of protein folding is a complex transition from a disordered to an ordered state. Here, we simulate a specific fast-folding protein at the point at which the native and denatured states are at equilibrium and show that obtaining an accurate description of the mechanisms of folding and unfolding is far from trivial. Using simple quantities which quantify the degree of native order is, in the case of this protein, clearly misleading. We show that an unbiased representation of the free-energy surface can be obtained; using such a representation we are able to redesign the landscape and thus modify, upon site-specific “mutations”, the folding and unfolding rates. This leads us to formulate a hypothesis to explain the very fast folding of many proteins.