Molecular dynamics simulations and membrane protein structure quality

Despite a growing repertoire of membrane protein structures (currently ∼120 unique structures), considerations of low resolution and crystallization in the absence of a lipid bilayer require the development of techniques to assess the global quality of membrane protein folds. This is also the case for assessment of, e.g. homology models of human membrane proteins based on structures of (distant) bacterial homologues. Molecular dynamics (MD) simulations may be used to help evaluate the quality of a membrane protein structure or model. We have used a structure of the bacterial ABC transporter MsbA which has the correct transmembrane helices but an incorrect handedness and topology of their packing to test simulation methods of quality assessment. An MD simulation of the MsbA model in a lipid bilayer is compared to a simulation of another bacterial ABC transporter, BtuCD. The latter structure has demonstrated good conformational stability in the same bilayer environment and over the same timescale (20 ns) as for the MsbA model simulation. A number of comparative analyses of the two simulations were performed to assess changes in the structural integrity of each protein. The results show a significant difference between the two simulations, chiefly due to the dramatic structural deformations of MsbA. We therefore propose that MD could become a useful quality control tool for membrane protein structural biology. In particular, it provides a way in which to explore the global conformational stability of a model membrane protein fold.     

Molecular dynamics simulations of membrane proteins

Membrane proteins control the traffic across cell membranes and thereby play an essential role in cell function from transport of various solutes to immune response via molecular recognition. Because it is very difficult to determine the structures of membrane proteins experimentally, computational methods have been increasingly used to study their structure and function. Here we focus on two classes of membrane proteins—ion channels and transporters—which are responsible for the generation of action potentials in nerves, muscles, and other excitable cells. We describe how computational methods have been used to construct models for these proteins and to study the transport mechanism. The main computational tool is the molecular dynamics (MD) simulation, which can be used for everything from refinement of protein structures to free energy calculations of transport processes. We illustrate with specific examples from gramicidin and potassium channels and aspartate transporters how the function of these membrane proteins can be investigated using MD simulations.     

Molecular dynamics simulations of the photoactive protein nitrile hydratase

Nitrile hydratase (NHase) is an enzyme used in the industrial biotechnological production of acrylamide. The active site, which contains nonheme iron or noncorrin cobalt, is buried in the protein core at the interface of two domains, α and β. Hydrogen bonds between βArg-56 and αCys-114 sulfenic acid (αCEA114) are important to maintain the enzymatic activity. The enzyme may be inactivated by endogenous nitric oxide (NO) and activated by absorption of photons of wavelength λ < 630 nm. To explain the photosensitivity and to propose structural determinants of catalytic activity, differences in the dynamics of light-active and dark-inactive forms of NHase were investigated using molecular dynamics (MD) modeling. To this end, a new set of force field parameters for nonstandard NHase active sites have been developed. The dynamics of the photodissociated NO ligand in the enzyme channel was analyzed using the locally enhanced sampling method, as implemented in the MOIL MD package. A series of 1 ns trajectories of NHases shows that the protonation state of the active site affects the dynamics of the catalytic water and NO ligand close to the metal center. MD simulations support the catalytic mechanism in which a water molecule bound to the metal ion directly attacks the nitrile carbon.     

Molecular dynamics simulations of the whey protein beta-lactoglobulin.

Molecular dynamics simulations have been used to model the motions and conformational behavior of the whey protein bovine beta-lactoglobulin. Simulations were performed for the protein by itself and complexed to a single retinol ligand located in a putative interior binding pocket. In the absence of the retinol ligand, the backbone loops around the opening of this interior pocket shifted inward to partially close off this cavity, similar to the shifts observed in previously reported molecular dynamics simulations of the uncomplexed form of the homologous retinol binding protein. The protein complexed with retinol does not exhibit the same conformational shifts. Conformational changes of this type could serve as a recognition signal allowing in vivo discrimination between the free and retinol complexed forms of the beta-lactoglobulin molecule. The unusual bending of the single alpha-helix observed in the simulations of retinol binding protein were not observed in the present calculations.     

Molecular dynamics simulations

Molecular dynamics simulations have become a standard tool for the investigation of biomolecules. Simulations are performed of ever bigger systems using more realistic boundary conditions and better sampling due to longer sampling times. Recently, realistic simulations of systems as complex as transmembrane channels have become feasible. Simulations aid our understanding of biochemical processes and give a dynamic dimension to structural data; for example, the transformation of harmless prion protein into the disease-causing agent has been modeled.     

Molecular dynamics simulations of membrane proteins under asymmetric ionic concentrations

A computational method is developed to allow molecular dynamics simulations of biomembrane systems under realistic ionic gradients and asymmetric salt concentrations while maintaining the conventional periodic boundary conditions required to minimize finite-size effects in an all-atom explicit solvent representation. The method, which consists of introducing a nonperiodic energy step acting on the ionic species at the edge of the simulation cell, is first tested with illustrative applications to a simple membrane slab model and a phospholipid membrane bilayer. The nonperiodic energy-step method is then used to calculate the reversal potential of the bacterial porin OmpF, a large cation-specific β-barrel channel, by simulating the I-V curve under an asymmetric 10:1 KCl concentration gradient. The calculated reversal potential of 28.6 mV is found to be in excellent agreement with the values of 26–27 mV measured from lipid bilayer experiments, thereby demonstrating that the method allows realistic simulations of nonequilibrium membrane transport with quantitative accuracy. As a final example, the pore domain of Kv1.2, a highly selective voltage-activated K⁺ channel, is simulated in a lipid bilayer under conditions that recreate, for the first time, the physiological K⁺ and Na⁺ concentration gradients and the electrostatic potential difference of living cells.     

Computer simulations of membrane protein folding: structure and dynamics

A lattice model of membrane proteins with a composite energy function is proposed to study their folding dynamics and native structures using Monte Carlo simulations. This model successfully predicts the seven helix bundle structure of sensory rhodopsin I by practicing a three-stage folding. Folding dynamics of a transmembrane segment into a helix is further investigated by varying the cooperativity in the formation of alpha helices for both random folding and assisted folding. The chain length dependence of the folding time of a hydrophobic segment to a helical state is studied for both free and anchored chains. An unusual length dependence in the folding time of anchored chains is observed.     

Molecular dynamics simulations of metalloproteins

Molecular dynamics simulations are now commonly applied to metalloproteins, despite the challenges introduced by the presence of metal ions. Force field parameters are nowadays available also for these 'exotic' atoms and several biological systems have been successfully studied. Some of the most relevant results and methodological advancements are reviewed.     

Molecular dynamics simulations investigation of neocarzinostatin chromophore-releasing pathways from the holo-NCS protein

The enediyne ring chromophore with strong DNA cleavage activity of neocarzinostatin is labile and therefore stabilization by forming the complex (carrying protein + chromophore: holo-NCS). Holo-NCS has gained much attention in clinical use as well as for drug delivery systems, but the chromophore-releasing mechanism to trigger binding to the target DNA with high affinity and producing DNA damage remain unclear. Three possible pathways were initially determined by conventional MD, essential dynamics and essential dynamics sampling. One of the paths runs along the naphthoate moiety; another runs along the amino sugar moiety; the third along the enediyne ring. Further, calculated forces and time by FPMD (force-probe molecular dynamics) suggest that the opening of the naphthoate moiety is most favorable pathway and Leu45, Phe76 and Phe78 all are key residues for chromophore release. In addition, conformational analyses indicate that the chromophore release is only local motions for the protein.     

Molecular dynamics simulations of xDNA

xDNA is a modified DNA, which contains natural as well as expanded bases. Expanded bases are generated by the addition of a benzene spacer to the natural bases. A set of AMBER force‐field parameters were derived for the expanded bases and the structural dynamics of the xDNA decamer ( xT5 ′ G xT A xC xG C xA xG T3′ ) · ( xA5′ C T xG C G xT A xC A3′) was explored using a 22 ns molecular dynamics simulation in explicit solvent. During the simulation, the duplex retained its Watson‐Crick base‐pairing and double helical structure, with deviations from the starting B‐form geometry towards A‐form; the deviations are mainly in the backbone torsion angles and in the helical parameters. The sugar pucker of the residues were distributed among a variety of modes; C2′ endo, C1′ exo, O4′ endo, C4′ exo, C2′ exo, and C3′ endo. The enhanced stacking interactions on account of the modification in the bases could help to retain the duplex nature of the helix with minor deviations from the ideal geometry. In our simulation, the xDNA showed a reduced minor groove width and an enlarged major groove width in comparison with the NMR structure. Both the grooves are larger than that of standard B‐DNA, but major groove width is larger than that of A‐DNA with almost equal minor groove width. The enlarged groove widths and the possibility of additional hydration in the grooves makes xDNA a potential molecule for various applications. © 2009 Wiley Periodicals, Inc. Biopolymers 91: 351–360, 2009. This article was originally published online as an accepted preprint. The “Published Online” date corresponds to the preprint version. You can request a copy of the preprint by emailing the Biopolymers editorial office at biopolymers@wiley.com     

Molecular dynamics simulations of carrabiose

Molecular mechanics calculations have been performed for the disaccharide carrabiose, one of the repeat units of β-carrageenan, as a general model for the (1→4)-linkage in the carrageenans. An adiabatic conformational energy map for this unsulfated molecule was prepared by constrained energy minimization and compared to a previously reported rigid-residue energy map for the sulfated molecule and to a similar adiabatic map for neocarrabiose, the related (1→3)-linked dimer repeat unit of β-carrageenan. Molecular dynamics simulations of this molecule in vacuo and in an aqueous (TIP3P) solution were calculated, and the observed motions were found to be generally consistent with the vacuum adiabatic energy map. Unlike the case observed in previous simulations of neocarrabiose, little salvation shift in the molecular conformation was observed for carrabiose. From the dynamics, the linkage was observed to be relatively flexible, as has been inferred from experiment on sulfated carrageenan polymers. © 1997 John Wiley & Sons, Inc.     

Molecular dynamics simulations as a tool for improving protein stability

Haloalkane dehalogenase (DhlA) was used as a model protein to explore the possibility to use molecular dynamics (MD) simulations as a tool to identify flexible regions in proteins that can serve as a target for stability enhancement by introduction of a disulfide bond. DhlA consists of two domains: an alpha/beta-hydrolase fold main domain and a cap domain composed of five alpha-helices. MD simulations of DhlA showed high mobility in a helix-loop-helix region in the cap domain, involving residues 184-211. A disulfide cross-link was engineered between residue 201 of this flexible region and residue 16 of the main domain. The mutant enzyme showed substantial changes in both thermal and urea denaturation. The oxidized form of the mutant enzyme showed an increase of the apparent transition temperature from 47.5 to 52.5 degrees C, whereas the T(m,app) of the reduced mutant decreased by more than 8 degrees C compared to the wild-type enzyme. Urea denaturation results showed a similar trend. Measurement of the kinetic stability showed that the introduction of the disulfide bond caused a decrease in activation free energy of unfolding of 0.43 kcal mol(-1) compared to the wild-type enzyme and also indicated that the helix-loop-helix region was involved early in the unfolding process. The results show that MD simulations are capable of identifying mobile protein domains that can successfully be used as a target for stability enhancement by the introduction of a disulfide cross-link.     

Molecular dynamics simulations of N-terminal peptides from a nucleotide binding protein

Molecular dynamics (MD) simulations of N-terminal peptides from lactate dehydrogenase (LDH) with increasing length and individual secondary structure elements were used to study their stability in relation to folding. Ten simulations of 1–2 ns of different peptides in water starting from the coordinates of the crystal structure were performed. The stability of the peptides was compared qualitatively by analyzing the root mean square deviation (RMSD) from the crystal structure, radius of gyration, secondary and tertiary structure, and solvent accessible surface area. In agreement with earlier MD studies, relatively short (< 15 amino acids) peptides containing individual secondary structure elements were generally found to be unstable; the hydrophobic α₁-helix of the nucleotide binding fold displayed a significantly higher stability, however. Our simulations further showed that the first βαβ supersecondary unit of the characteristic dinucleotide binding fold (Rossmann fold) of LDH is somewhat more stable than other units of similar length and that the α₂-helix, which unfolds by itself, is stabilized by binding to this unit. This finding suggests that the first βαβ unit could function as an N-terminal folding nucleus, upon which the remainder of the polypeptide chain can be assembled. Indeed, simulations with longer units (βαβα and βαβαββ) showed that all structural elements of these units are rather stable. The outcome of our studies is in line with suggestions that folding of the N-terminal portion of LDH in vivo can be a cotranslational process that takes place during the ribosomal peptide synthesis.