 |
|
 |
 | We have recently changed our website. Looking for the previous site? |
Abstracts of Articles in 2007
To request any journal reprints, please contact us.- Flap opening dynamics in HIV-1 protease explored with a coarse-grained model
- Optimizing the Poisson Dielectric Boundary with Explicit Solvent Forces and Energies: Lessons Learned with Atom-Centered Dielectric Functions
- Generalized Born model with a simple, robust molecular volume correction
- Unraveling the Three-Metal-Ion Catalytic Mechanism of the DNA Repair Enzyme Endonuclease IV
- Towards an in vivo biologically inspired nanofactory
- SR protein kinase 1 is resilient to inactivation
- Binding Pathways of Ligands to HIV-1 Protease: Coarse-grained and Atomistic Simulations
- HIV-1 protease substrate binding and product release pathways explored with coarse-grained molecular dynamics
- Structure and Dynamics of the Full-Length Lipid-Modified H-Ras Protein in a 1,2-dimyristoylglycero-3-phosphocholine Bilayer
- Continuum simulations of acetylcholine diffusion with reaction-determined boundaries in neuromuscular junction models
- Finite Element Analysis of the Time-Dependent Smoluchowski Equation for Acetylcholinesterase Reaction Rate Calculations
- Peptide insertion, positioning, and stabilization in a membrane: Insight from an all-atom molecular dynamics simulation
- Comparative MD analysis of the stability of Transthyretin providing insight into the fibrillation mechanism
- Multivariate analysis of conserved sequence-structure relationships in kinesins: Coupling of the active site and a tubulin-binding domain
- Molecular Dynamics Simulations of Metalloproteinases types 2 and 3 Reveal Differences in the Dynamic Behavior of the S1' Binding Pocket
- Improved Boundary Element Methods for Poisson-Boltzmann Electrostatic Potential and Force Calculations
- Acetylcholinesterase: Mechanisms of Covalent Inhibition of Wild-Type and H447I Mutant Determined by Computational Analyses
- Barriers to Ion Translocation in Cationic and Anionic Receptors from the Cys-Loop Family
- Dealing with bound waters in a site: Do they leave or stay?
- Remarkable Loop Flexibility in Avian Influenza N1 and its Implications for Antiviral Drug Design
- “New-version-fast-multipole-method” accelerated electrostatic calculations in biomolecular systems
- Application of the level-set method to the implicit solvation of nonpolar molecules
- Nanosecond-Timescale Conformational Dynamics of the Human α7 Nicotinic Acetylcholine Receptor
- Functional and Structural Insights Revealed by Molecular Dynamics Simulations of an Essential RNA Editing Ligase in Trypanosoma brucei
- Free Energy Profile of H-ras Membrane Anchor upon Membrane Insertion
- Generalized gradient-augmented harmonic Fourier beads method with multiple atomic and/or center-of-mass positional restraints
- Solvent reaction field potential inside an uncharged globular protein: A bridge between implicit and explicit solvent models?
- H-ras Protein in a Bilayer: Interaction and Structure Perturbation
- A proposed signaling motif for nuclear import in mRNA processing via the formation of Arginine Claw
- Electrodiffusion: A continuum modeling framework for biomolecular systems with realistic spatiotemporal resolution
- Dynamics, Hydration, and Motional Averaging of a Loop-Gated Artificial Protein Cavity: The W191G Mutant of Cytochrome c Peroxidase in Water as Revealed by Molecular Dynamics Simulations
- Sampling of slow diffusive conformational transitions with accelerated molecular dynamics
- Estimating Kinetic Rates from Accelerated Molecular Dynamics Simulations: Alanine Dipeptide in Explicit Solvent as a Case Study
- Accelerated entropy estimates with accelerated dynamics
Flap opening dynamics in HIV-1 protease explored with a coarse-grained modelJournal of Structural Biology, Vol. 157, Issue 3, pp. 606-615 (2007) [PubMed 17029846]
We present a one-bead coarse-grained model that enables dynamical simulations of proteins on the time scale of tens of microseconds. The parameterization of the force field includes accurate conformational terms that allow for fast and reliable exploration of the configurational space. The model is applied to the dynamics of flap opening in HIV-1 protease. The experimental structure of the recently crystallized semi-open conformation of HIV-1 protease is well reproduced in the simulation, which supports the accuracy of our model. Thanks to very long simulations and extensive sampling of opening and closing events, we also investigate the thermodynamics and kinetics of the opening process. We have shown that the effect of the solvent slows down the dynamics to the experimentally observed time scales. The model is found to be reliable for application to substrate docking simulations, which are currently in progress.
We present a one-bead coarse-grained model that enables dynamical simulations of proteins on the time scale of tens of microseconds. The parameterization of the force field includes accurate conformational terms that allow for fast and reliable exploration of the configurational space. The model is applied to the dynamics of flap opening in HIV-1 protease. The experimental structure of the recently crystallized semi-open conformation of HIV-1 protease is well reproduced in the simulation, which supports the accuracy of our model. Thanks to very long simulations and extensive sampling of opening and closing events, we also investigate the thermodynamics and kinetics of the opening process. We have shown that the effect of the solvent slows down the dynamics to the experimentally observed time scales. The model is found to be reliable for application to substrate docking simulations, which are currently in progress.
Optimizing the Poisson Dielectric Boundary with Explicit Solvent Forces and Energies: Lessons Learned with Atom-Centered Dielectric FunctionsJournal of Chemical Theory and Computation, Vol. 3, No. 1, pp. 170-183 (2007)
Accurate implicit solvent models require parameters that have been optimized in a system- and/or atom-specific manner based on experimental data or more rigorous explicit solvent simulations. Models based on the Poisson or Poisson-Boltzmann equation are particularly sensitive to the nature and location of the boundary which separates the low dielectric solute from the high dielectric solvent. Here we present a novel method for optimizing the solute radii, which define the dielectric boundary, based on forces and energies from explicit solvent simulations. We use this method to optimize radii for protein systems defined by AMBER /ff/99 partial charges and a spline-smoothed solute surface. The spline-smoothed surface is an atom-centered dielectric function that enables stable and efficient force calculations. We explore the relative performance of radii optimized with forces alone and those optimized with forces and energies. We show that our radii reproduce the explicit solvent forces and energies more accurately than four other parameter sets commonly used in conjunction with the AMBER force field, each of which has been appropriately scaled for spline-smoothed surfaces. Finally, we demonstrate that spline-smoothed surfaces show surprising accuracy for small, compact systems, but may have limitations for highly-solvated protein systems. The optimization method presented here is efficient and applicable to any system with explicit solvent parameters. It can be used to determine the optimal continuum parameters when experimental solvation energies are unavailable and the computational costs of explicit solvent charging free energies are prohibitive.
Accurate implicit solvent models require parameters that have been optimized in a system- and/or atom-specific manner based on experimental data or more rigorous explicit solvent simulations. Models based on the Poisson or Poisson-Boltzmann equation are particularly sensitive to the nature and location of the boundary which separates the low dielectric solute from the high dielectric solvent. Here we present a novel method for optimizing the solute radii, which define the dielectric boundary, based on forces and energies from explicit solvent simulations. We use this method to optimize radii for protein systems defined by AMBER /ff/99 partial charges and a spline-smoothed solute surface. The spline-smoothed surface is an atom-centered dielectric function that enables stable and efficient force calculations. We explore the relative performance of radii optimized with forces alone and those optimized with forces and energies. We show that our radii reproduce the explicit solvent forces and energies more accurately than four other parameter sets commonly used in conjunction with the AMBER force field, each of which has been appropriately scaled for spline-smoothed surfaces. Finally, we demonstrate that spline-smoothed surfaces show surprising accuracy for small, compact systems, but may have limitations for highly-solvated protein systems. The optimization method presented here is efficient and applicable to any system with explicit solvent parameters. It can be used to determine the optimal continuum parameters when experimental solvation energies are unavailable and the computational costs of explicit solvent charging free energies are prohibitive.
Generalized Born model with a simple, robust molecular volume correctionJournal of Chemical Theory and Computation, Vol. 3, No. 1, pp. 156-169 (2007) [PubMed 21072141]
Generalized Born (GB) models provide a computationally efficient means of representing the electrostatic effects of solvent and are widely used, especially in molecular dynamics (MD). A class of particularly fast GB models is based on integration over an interior volume approximated as a pairwise union of atom spheres -- effectively, the interior is defined by a van der Waals rather than Lee-Richards molecular surface. The approximation is computationally efficient, but if uncorrected, allows for high dielectric (water) regions smaller than a water molecule between atoms, leading to decreased accuracy. Here, an earlier pairwise GB model is extended by a simple analytic correction term that largely alleviates the problem by correctly describing the solvent-excluded volume of each pair of atoms. The correction term introduces a free energy barrier to the separation of non-bonded atoms. This free energy barrier is seen in explicit solvent and Lee-Richards molecular surface implicit solvent calculations, but has been absent from earlier pairwise GB models. When used in MD, the correction term yields protein hydrogen bond length distributions and polypeptide conformational ensembles that are in better agreement with explicit solvent results than earlier pairwise models. The robustness and simplicity of the correction preserves the efficiency of the pairwise GB models while making them a better approximation to reality.
Generalized Born (GB) models provide a computationally efficient means of representing the electrostatic effects of solvent and are widely used, especially in molecular dynamics (MD). A class of particularly fast GB models is based on integration over an interior volume approximated as a pairwise union of atom spheres -- effectively, the interior is defined by a van der Waals rather than Lee-Richards molecular surface. The approximation is computationally efficient, but if uncorrected, allows for high dielectric (water) regions smaller than a water molecule between atoms, leading to decreased accuracy. Here, an earlier pairwise GB model is extended by a simple analytic correction term that largely alleviates the problem by correctly describing the solvent-excluded volume of each pair of atoms. The correction term introduces a free energy barrier to the separation of non-bonded atoms. This free energy barrier is seen in explicit solvent and Lee-Richards molecular surface implicit solvent calculations, but has been absent from earlier pairwise GB models. When used in MD, the correction term yields protein hydrogen bond length distributions and polypeptide conformational ensembles that are in better agreement with explicit solvent results than earlier pairwise models. The robustness and simplicity of the correction preserves the efficiency of the pairwise GB models while making them a better approximation to reality.
Unraveling the Three-Metal-Ion Catalytic Mechanism of the DNA Repair Enzyme Endonuclease IVProceedings of the National Academy of Sciences of the USA, Vol. 104, No. 5, pp. 1465-1470 (2007) [PubMed 17242363]
Endonuclease IV belongs to a class of important apurinic/apyrimidinic (AP) endonucleases involved in DNA repair. Although a structure-based mechanistic hypothesis has been put forth for this enzyme, the detailed catalytic mechanism has remained unknown. Using thermodynamic integration in the context of ab initio QM/MM molecular dynamics (AIMD), we examined certain aspects of the phosphodiester cleavage step in the mechanism. We found the reaction proceeded through a synchronous bimolecular (ANDN) mechanism with reaction free energy and barrier of -3.5 and 20.6 kcal/mol, in agreement with experimental estimates. In the course of the reaction the tri-nuclear active site of endonuclease IV underwent dramatic local conformational changes -- shifts in the mode of coordination of both substrate and first-shell ligands. This qualitative finding supports the notion that structural rearrangements in the active sites of multinuclear enzymes are integral to biological function.
Endonuclease IV belongs to a class of important apurinic/apyrimidinic (AP) endonucleases involved in DNA repair. Although a structure-based mechanistic hypothesis has been put forth for this enzyme, the detailed catalytic mechanism has remained unknown. Using thermodynamic integration in the context of ab initio QM/MM molecular dynamics (AIMD), we examined certain aspects of the phosphodiester cleavage step in the mechanism. We found the reaction proceeded through a synchronous bimolecular (ANDN) mechanism with reaction free energy and barrier of -3.5 and 20.6 kcal/mol, in agreement with experimental estimates. In the course of the reaction the tri-nuclear active site of endonuclease IV underwent dramatic local conformational changes -- shifts in the mode of coordination of both substrate and first-shell ligands. This qualitative finding supports the notion that structural rearrangements in the active sites of multinuclear enzymes are integral to biological function.
Towards an in vivo biologically inspired nanofactoryNature Nanotechnology, Vol. 2, No. 1, pp. 3-7 (2007) [PubMed 18654192]
Nanotechnology is having a major impact on medicine and the treatment of disease, notably in imaging and targeted drug delivery. It may, however, be possible to go even further and design 'pseudo-cell' nanofactories that work with molecules already in the body to fight disease.
Nanotechnology is having a major impact on medicine and the treatment of disease, notably in imaging and targeted drug delivery. It may, however, be possible to go even further and design 'pseudo-cell' nanofactories that work with molecules already in the body to fight disease.
SR protein kinase 1 is resilient to inactivationStructure, Vol. 15, Issue 1, pp. 123-133 (2007) [PubMed 17223538]
SR protein kinase 1 (SRPK1) is a constitutively active kinase, which processively phosphorylates multiple serines within its substrates, ASF/SF2. We describe crystallographic, molecular dynamics, and biochemical results that shed light on how SRPK1 preserves its constitutive active conformation. Our structure reveals that unlike other known active kinase structures, the activation loop remains in an active state without any specific intraprotein interactions. Moreover, SRPK1 remains active despite extensive mutation to the activation segment. Molecular dynamics simulations reveal that SRPK1 partially absorbs the effect of mutations by forming compensatory interactions that maintain a catalytically competent chemical environment. Furthermore, SRPK1 is similarly resistant to deletion of its spacer loop region. Based upon a model of SRPK1 bound to a segment encompassing the docking motif and active-site peptide of ASF/SF2, we suggest a mechanism for processive phosphorylation and propose that the atypical resiliency we observed is critical for SRPK1's processive activity.
SR protein kinase 1 (SRPK1) is a constitutively active kinase, which processively phosphorylates multiple serines within its substrates, ASF/SF2. We describe crystallographic, molecular dynamics, and biochemical results that shed light on how SRPK1 preserves its constitutive active conformation. Our structure reveals that unlike other known active kinase structures, the activation loop remains in an active state without any specific intraprotein interactions. Moreover, SRPK1 remains active despite extensive mutation to the activation segment. Molecular dynamics simulations reveal that SRPK1 partially absorbs the effect of mutations by forming compensatory interactions that maintain a catalytically competent chemical environment. Furthermore, SRPK1 is similarly resistant to deletion of its spacer loop region. Based upon a model of SRPK1 bound to a segment encompassing the docking motif and active-site peptide of ASF/SF2, we suggest a mechanism for processive phosphorylation and propose that the atypical resiliency we observed is critical for SRPK1's processive activity.
Binding Pathways of Ligands to HIV-1 Protease: Coarse-grained and Atomistic SimulationsChemical Biology & Drug Design, Vol. 69, Issue 1, pp. 5-13 (2007) [PubMed 17313452]
Multiscale simulations (coarse-grained Brownian dynamics simulations and all-atom molecular dynamics simulations in implicit solvent) were applied to reveal the binding processes of ligands as they enter the binding site of the HIV-1 protease. The initial structures used for the molecular dynamics simulations were generated based on the Brownian dynamics trajectories, and this is the first molecular dynamics simulation of modeling the association of a ligand with the protease. We found that a protease substrate successfully binds to the protein when the flaps are fully open. Surprisingly, a smaller cyclic urea inhibitor (XK263) can reach the binding site when the flaps are not fully open. However, if the flaps are nearly closed, the inhibitor must rearrange or binding can fail because the inhibitor cannot attain proper conformations to enter the binding site. Both the peptide substrate and XK263 can also affect the protein's internal motion, which may help the flaps to open. Simulations allows us to efficiently study the ligand binding processes and may help those who study drug discovery to find optimal association pathways and to design those ligands with the best binding kinetics.
Multiscale simulations (coarse-grained Brownian dynamics simulations and all-atom molecular dynamics simulations in implicit solvent) were applied to reveal the binding processes of ligands as they enter the binding site of the HIV-1 protease. The initial structures used for the molecular dynamics simulations were generated based on the Brownian dynamics trajectories, and this is the first molecular dynamics simulation of modeling the association of a ligand with the protease. We found that a protease substrate successfully binds to the protein when the flaps are fully open. Surprisingly, a smaller cyclic urea inhibitor (XK263) can reach the binding site when the flaps are not fully open. However, if the flaps are nearly closed, the inhibitor must rearrange or binding can fail because the inhibitor cannot attain proper conformations to enter the binding site. Both the peptide substrate and XK263 can also affect the protein's internal motion, which may help the flaps to open. Simulations allows us to efficiently study the ligand binding processes and may help those who study drug discovery to find optimal association pathways and to design those ligands with the best binding kinetics.
HIV-1 protease substrate binding and product release pathways explored with coarse-grained molecular dynamicsBiophysical Journal, Vol. 92, Issue 12, pp. 4179-4187 (2007) [PubMed 17384072]
We analyze the encounter of a peptide substrate with the native HIV-1 protease, the mechanism of substrate incorporation in the binding cleft, and the dissociation of products after substrate hydrolysis. To account for the substrate, we extend a coarse-grained model force field which we previously developed to study the flap opening dynamics of HIV-1 protease on a microsecond time scale. Molecular and Langevin dynamics simulations show that the flaps need to open for the peptide to bind and that the protease interaction with the substrate influences the flap opening frequency and interval. On the other hand, release of the products does not require flap opening because they can slide out from the binding cleft to the sides of the enzyme. Our data show that in the protease-substrate complex the highest fluctuations correspond to the 17- and 39-turns and the substrate motion is anti-correlated with the 39-turn. Moreover, the active site residues and the flap tips move in phase with the peptide. We suggest some mechanistic principles for how the flexibility of the protein may be involved in ligand binding and release.
We analyze the encounter of a peptide substrate with the native HIV-1 protease, the mechanism of substrate incorporation in the binding cleft, and the dissociation of products after substrate hydrolysis. To account for the substrate, we extend a coarse-grained model force field which we previously developed to study the flap opening dynamics of HIV-1 protease on a microsecond time scale. Molecular and Langevin dynamics simulations show that the flaps need to open for the peptide to bind and that the protease interaction with the substrate influences the flap opening frequency and interval. On the other hand, release of the products does not require flap opening because they can slide out from the binding cleft to the sides of the enzyme. Our data show that in the protease-substrate complex the highest fluctuations correspond to the 17- and 39-turns and the substrate motion is anti-correlated with the 39-turn. Moreover, the active site residues and the flap tips move in phase with the peptide. We suggest some mechanistic principles for how the flexibility of the protein may be involved in ligand binding and release.
Structure and Dynamics of the Full-Length Lipid-Modified H-Ras Protein in a 1,2-dimyristoylglycero-3-phosphocholine BilayerJournal of Medicinal Chemistry, Vol. 50, No. 4, pp. 674-684 (2007) [PubMed 17263520]
Ras proteins regulate signal transduction processes that control cell growth and proliferation. Their dysregulation is a common cause of human tumors. Atomic level structural and dynamical information in a membrane environment is crucial for understanding signaling specificity among Ras isoforms and for the design of selective anti-cancer agents. Here, the structure of the full-length H-Ras protein in complex with a 1,2-dimyristoylglycero- 3-phosphocholine (DMPC) bilayer obtained from modeling and all-atom explicit solvent molecular dynamics simulations, as well as experimental validation of the main results, are presented. We find that, in addition to the lipid anchor, H-Ras membrane binding involves direct interaction of residues in the catalytic domain with DMPC phosphates. Two modes of binding (possibly modulated by GTP/GDP exchange) differing in the orientation and bilayer contact of the soluble domain as well as in the participation of the flexible linker in membrane binding are proposed. These results are supported by our initial in vivo experiments. The overall structure of the protein and the bilayer remain similar to that of the isolated components, with few localized structural and dynamical changes. The implications of the results to membrane lateral segregation and other aspects of Ras signaling are discussed.
Ras proteins regulate signal transduction processes that control cell growth and proliferation. Their dysregulation is a common cause of human tumors. Atomic level structural and dynamical information in a membrane environment is crucial for understanding signaling specificity among Ras isoforms and for the design of selective anti-cancer agents. Here, the structure of the full-length H-Ras protein in complex with a 1,2-dimyristoylglycero- 3-phosphocholine (DMPC) bilayer obtained from modeling and all-atom explicit solvent molecular dynamics simulations, as well as experimental validation of the main results, are presented. We find that, in addition to the lipid anchor, H-Ras membrane binding involves direct interaction of residues in the catalytic domain with DMPC phosphates. Two modes of binding (possibly modulated by GTP/GDP exchange) differing in the orientation and bilayer contact of the soluble domain as well as in the participation of the flexible linker in membrane binding are proposed. These results are supported by our initial in vivo experiments. The overall structure of the protein and the bilayer remain similar to that of the isolated components, with few localized structural and dynamical changes. The implications of the results to membrane lateral segregation and other aspects of Ras signaling are discussed.
Continuum simulations of acetylcholine diffusion with reaction-determined boundaries in neuromuscular junction modelsBiophysical Chemistry, Vol. 127, Issue 3, pp. 129-139 (2007) [PubMed 17307283]
The reaction-diffusion system of the neuromuscular junction has been modeled in 3D using the finite element package FEtk. The numerical solution of the dynamics of acetylcholine with the detailed reaction processes of acetylcholinesterases and nicotinic acetylcholine receptors has been discussed with the reaction-determined boundary conditions. The simulation results describe the detailed acetylcholine hydrolysis process, and reveal the timedependent interconversion of the closed and open states of the acetylcholine receptors as well as the percentages of unliganded/monoliganded/diliganded states during the neurotransmission. The finite element method has demonstrated its flexibility and robustness in modeling large biological systems.
The reaction-diffusion system of the neuromuscular junction has been modeled in 3D using the finite element package FEtk. The numerical solution of the dynamics of acetylcholine with the detailed reaction processes of acetylcholinesterases and nicotinic acetylcholine receptors has been discussed with the reaction-determined boundary conditions. The simulation results describe the detailed acetylcholine hydrolysis process, and reveal the timedependent interconversion of the closed and open states of the acetylcholine receptors as well as the percentages of unliganded/monoliganded/diliganded states during the neurotransmission. The finite element method has demonstrated its flexibility and robustness in modeling large biological systems.
Finite Element Analysis of the Time-Dependent Smoluchowski Equation for Acetylcholinesterase Reaction Rate CalculationsBiophysical Journal, Vol. 92, Issue 10, pp. 3397-3406 (2007) [PubMed 17307827]
This article describes the numerical solution of the time-dependent Smoluchowski equation to study diffusion in biomolecular systems. Specifically, finite element methods have been developed to calculate ligand binding rate constants for large biomolecules. The resulting software has been validated and applied to the mouse acetylcholinesterase monomer and several tetramers. Rates for inhibitor binding to mAChE were calculated at various ionic strengths with several different time steps. Calculated rates show very good agreement with experimental and theoretical steady-state studies. Furthermore, these finite element methods require significantly fewer computational resources than existing particle-based Brownian dynamics methods and are robust for complicated geometries. The key finding of biological importance is that the rate accelerations of the monomeric and tetrameric mAChE that result from electrostatic steering are preserved under the non-steadystate conditions that are expected to occur in physiological circumstances.
This article describes the numerical solution of the time-dependent Smoluchowski equation to study diffusion in biomolecular systems. Specifically, finite element methods have been developed to calculate ligand binding rate constants for large biomolecules. The resulting software has been validated and applied to the mouse acetylcholinesterase monomer and several tetramers. Rates for inhibitor binding to mAChE were calculated at various ionic strengths with several different time steps. Calculated rates show very good agreement with experimental and theoretical steady-state studies. Furthermore, these finite element methods require significantly fewer computational resources than existing particle-based Brownian dynamics methods and are robust for complicated geometries. The key finding of biological importance is that the rate accelerations of the monomeric and tetrameric mAChE that result from electrostatic steering are preserved under the non-steadystate conditions that are expected to occur in physiological circumstances.
Peptide insertion, positioning, and stabilization in a membrane: Insight from an all-atom molecular dynamics simulationBiopolymers, Vol. 85, Issue 5-6, pp. 490-497 (2007)
Peptide insertion, positioning, and stabilization in a model membrane are probed via an all-atom molecular dynamics simulation. One peptide (WL5) is simulated in each leaflet of a solvated dimyristoylglycero-3-phosphate (DMPC) membrane. Within the first 5 ns, the peptides spontaneously insert into the membrane and then stabilize during the remaining 70 ns of simulation time. In both leaflets, the peptides localize to the membrane interface, and this localization is attributed to the formation of peptide-lipid hydrogen bonds. We show that the single tryptophan residue in each peptide contributes significantly to these hydrogen bonds; specifically, the nitrogen heteroatom of the indole ring plays a critical role. The tilt angles of the indole rings relative to the membrane normal in the upper and lower leaflets are approximately 26° and 54°, respectively. The tilt angles of the entire peptide chain are 62° and 74°. The membrane induces conformations of the peptide that are characteristic of β-sheets, and the peptide enhances the lipid ordering in the membrane. Finally, the diffusion rate of the peptides in the membrane plane is calculated (based on experimental peptide concentrations) to be approximately 6 Å2/ns, thus suggesting a 500 ns time scale for intermolecular interactions.
Peptide insertion, positioning, and stabilization in a model membrane are probed via an all-atom molecular dynamics simulation. One peptide (WL5) is simulated in each leaflet of a solvated dimyristoylglycero-3-phosphate (DMPC) membrane. Within the first 5 ns, the peptides spontaneously insert into the membrane and then stabilize during the remaining 70 ns of simulation time. In both leaflets, the peptides localize to the membrane interface, and this localization is attributed to the formation of peptide-lipid hydrogen bonds. We show that the single tryptophan residue in each peptide contributes significantly to these hydrogen bonds; specifically, the nitrogen heteroatom of the indole ring plays a critical role. The tilt angles of the indole rings relative to the membrane normal in the upper and lower leaflets are approximately 26° and 54°, respectively. The tilt angles of the entire peptide chain are 62° and 74°. The membrane induces conformations of the peptide that are characteristic of β-sheets, and the peptide enhances the lipid ordering in the membrane. Finally, the diffusion rate of the peptides in the membrane plane is calculated (based on experimental peptide concentrations) to be approximately 6 Å2/ns, thus suggesting a 500 ns time scale for intermolecular interactions.
Comparative MD analysis of the stability of Transthyretin providing insight into the fibrillation mechanismBiopolymers, Vol. 86, Issue 1, pp. 73-82 (2007) [PubMed 17315201]
Proteins can misfold and aggregate, which is believed to be the cause of a variety of diseases, affecting very diverse organs in the body. Many questions about the nature of aggregation and the proteins that are involved in these events are still left unanswered. One of the proteins that is known to form amyloids is transthyretin (TTR), the secondary transporter of thyroxine, and transporter of retinol-binding protein. Several experimental results have helped to explain this aberrant behavior of TTR; however, structural insights of the amyloidgenic process are still lacking. Therefore, we have used all-atom MD simulation and free energy calculations to study the initial phase of this process. We have calculated the free energy changes of the initial tetramer dissociation under different conditions and in the presence of thyroxine. We show that tetramer formation is indeed only thermodynamically favorable in neutral pH conditions. We find that binding of two thyroxine molecules stabilizes the complex, and that this occurs with negative cooperativity. In addition to the energetic calculations, we have also investigated the dominant motions of the TTR and found that only the dimeric form of the protein could undergo the initial fibril formation.
Proteins can misfold and aggregate, which is believed to be the cause of a variety of diseases, affecting very diverse organs in the body. Many questions about the nature of aggregation and the proteins that are involved in these events are still left unanswered. One of the proteins that is known to form amyloids is transthyretin (TTR), the secondary transporter of thyroxine, and transporter of retinol-binding protein. Several experimental results have helped to explain this aberrant behavior of TTR; however, structural insights of the amyloidgenic process are still lacking. Therefore, we have used all-atom MD simulation and free energy calculations to study the initial phase of this process. We have calculated the free energy changes of the initial tetramer dissociation under different conditions and in the presence of thyroxine. We show that tetramer formation is indeed only thermodynamically favorable in neutral pH conditions. We find that binding of two thyroxine molecules stabilizes the complex, and that this occurs with negative cooperativity. In addition to the energetic calculations, we have also investigated the dominant motions of the TTR and found that only the dimeric form of the protein could undergo the initial fibril formation.
Multivariate analysis of conserved sequence-structure relationships in kinesins: Coupling of the active site and a tubulin-binding domainJournal of Molecular Biology, Vol. 368, Issue 5, pp. 1231-1248 (2007) [PubMed 17399740]
An extensive computational analysis of available sequence and crystal structure data was used to identify functionally important residue interactions within the motor domain of the kinesin molecular motor. Principal component analysis revealed that all current kinesin crystal structures reside in one of two main conformations, which differ at the active site, and in the position of a microtubule-binding sub-domain relative to a rigid central core. This sub-domain consists of secondary structure elements α4-loop12-α5-loop13 and contains a conserved hydrophilic surface patch that may be involved in strong binding to microtubules. A hinge point for the sub-domain motion lies near a conserved glycine at position 292. Statistical coupling analysis revealed a network of co-evolving positions that link this region to the nucleotide-binding site, via a highly conserved histidine in the switch I loop. The data are consistent with a model in which the nucleotide status of the active site shifts kinesin between weak and strong binding conformations via reconfiguration of the identified sub-domain. Our data provide a statistically supported framework for further examination of this and other structure-function relationships in the kinesin family.
An extensive computational analysis of available sequence and crystal structure data was used to identify functionally important residue interactions within the motor domain of the kinesin molecular motor. Principal component analysis revealed that all current kinesin crystal structures reside in one of two main conformations, which differ at the active site, and in the position of a microtubule-binding sub-domain relative to a rigid central core. This sub-domain consists of secondary structure elements α4-loop12-α5-loop13 and contains a conserved hydrophilic surface patch that may be involved in strong binding to microtubules. A hinge point for the sub-domain motion lies near a conserved glycine at position 292. Statistical coupling analysis revealed a network of co-evolving positions that link this region to the nucleotide-binding site, via a highly conserved histidine in the switch I loop. The data are consistent with a model in which the nucleotide status of the active site shifts kinesin between weak and strong binding conformations via reconfiguration of the identified sub-domain. Our data provide a statistically supported framework for further examination of this and other structure-function relationships in the kinesin family.
Molecular Dynamics Simulations of Metalloproteinases types 2 and 3 Reveal Differences in the Dynamic Behavior of the S1' Binding PocketCurrent Pharmaceutical Design, Vol. 13, No. 34, pp. 3471-3475 (2007) [PubMed 18220784]
Matrix Metalloproteinases (MMPs) are zinc-containing proteinases that are responsible for the metabolism of extracellular matrix proteins. Overexpression of MMPs has been associated with a wide range of pathological diseases such as arthritis, cancer, multiple sclerosis and Alzheimer's disease. The excessive and unregulated activity of Matrix Metalloproteinases type 2 (MMP-2), also known as gelantinase A, has been identified in a number of cancer metastases. Several MMP inhibitors (MMPi) have been proposed in the literature aiming to interfere in the MMPs activity. In this work we performed long MD simulations in order to study the dynamical behavior of the binding pocket S1' in the apo forms of MMP type 2 and 3, and identify, at the molecular level, the structural properties relevant for the designing of specific inhibitor MMP-2.
Matrix Metalloproteinases (MMPs) are zinc-containing proteinases that are responsible for the metabolism of extracellular matrix proteins. Overexpression of MMPs has been associated with a wide range of pathological diseases such as arthritis, cancer, multiple sclerosis and Alzheimer's disease. The excessive and unregulated activity of Matrix Metalloproteinases type 2 (MMP-2), also known as gelantinase A, has been identified in a number of cancer metastases. Several MMP inhibitors (MMPi) have been proposed in the literature aiming to interfere in the MMPs activity. In this work we performed long MD simulations in order to study the dynamical behavior of the binding pocket S1' in the apo forms of MMP type 2 and 3, and identify, at the molecular level, the structural properties relevant for the designing of specific inhibitor MMP-2.
Improved Boundary Element Methods for Poisson-Boltzmann Electrostatic Potential and Force CalculationsJournal of Chemical Theory and Computation, Vol. 3, Issue 3, pp. 1134-1142 (2007)
A patch representation differing from the traditional treatments in the boundary element method (BEM) is presented, which we call the constant "node patch" method. Its application to solving the Poisson-Boltzmann equation (PBE) demonstrates considerable improvement in speed compared with the constant element and linear element methods. In addition, for the node-based BEMs, we propose an efficient interpolation method for the calculation of the electrostatic stress tensor and PB force on the solvated molecular surface. This force calculation is simply an O(N) algorithm (N is the number of elements). Moreover, our calculations also show that the geometric factor correction in the boundary integral equations (BIEs) significantly increases the accuracy of the potential solution on the boundary, and thereby the PB force calculation.
A patch representation differing from the traditional treatments in the boundary element method (BEM) is presented, which we call the constant "node patch" method. Its application to solving the Poisson-Boltzmann equation (PBE) demonstrates considerable improvement in speed compared with the constant element and linear element methods. In addition, for the node-based BEMs, we propose an efficient interpolation method for the calculation of the electrostatic stress tensor and PB force on the solvated molecular surface. This force calculation is simply an O(N) algorithm (N is the number of elements). Moreover, our calculations also show that the geometric factor correction in the boundary integral equations (BIEs) significantly increases the accuracy of the potential solution on the boundary, and thereby the PB force calculation.
Acetylcholinesterase: Mechanisms of Covalent Inhibition of Wild-Type and H447I Mutant Determined by Computational AnalysesJournal of the American Chemical Society, Vol. 129, Issue 20, pp. 6562-6570 (2007) [PubMed 17461584]
The reaction mechanisms of two inhibitors TFK+ and TFK0 binding to both the wild-type and H447I mutant mouse acetylcholinesterase (mAChE) have been investigated by using a combined ab initio quantum mechanical/molecular mechanical (QM/MM) approach and classical molecular dynamics (MD) simulations. In the wild-type mAChE, the binding reactions of TFK+ and TFK0 are both spontaneous processes, which proceed through the nucleophilic addition of the Ser203-Oγ to the carbonyl-C of TFK+ or TFK0, accompanied with a simultaneous proton transfer from Ser203 to His447. No barrier is found along the reaction paths, consistent with the experimental reaction rates approaching the diffusion-controlled limit. By contrast, TFK+ binding to the H447I mutant may proceed with a different reaction mechanism. A water molecule takes over the role of His447 and participates in the bond breaking and forming as a "charge relayer". Unlike in the wild-type mAChE case, Glu334, a conserved residue from the catalytic triad, acts as a catalytic base in the reaction. The calculated energy barrier for this reaction is about 8 kcal/mol. These predictions await experimental verification. In the case of the neutral ligand TFK0, however, multiple MD simulations on the TFK0/H447I complex reveal that none of the water molecules can be retained in the active site as a "catalytic" water. Furthermore, our alchemical free energy calculation also suggests that the binding of TFK0 to H447I is much weaker than that of TFK+. Taken together, our computational studies confirm that TFK0 is almost inactive in the H447I mutant and also provide detailed mechanistic insights into the experimental observations.
The reaction mechanisms of two inhibitors TFK+ and TFK0 binding to both the wild-type and H447I mutant mouse acetylcholinesterase (mAChE) have been investigated by using a combined ab initio quantum mechanical/molecular mechanical (QM/MM) approach and classical molecular dynamics (MD) simulations. In the wild-type mAChE, the binding reactions of TFK+ and TFK0 are both spontaneous processes, which proceed through the nucleophilic addition of the Ser203-Oγ to the carbonyl-C of TFK+ or TFK0, accompanied with a simultaneous proton transfer from Ser203 to His447. No barrier is found along the reaction paths, consistent with the experimental reaction rates approaching the diffusion-controlled limit. By contrast, TFK+ binding to the H447I mutant may proceed with a different reaction mechanism. A water molecule takes over the role of His447 and participates in the bond breaking and forming as a "charge relayer". Unlike in the wild-type mAChE case, Glu334, a conserved residue from the catalytic triad, acts as a catalytic base in the reaction. The calculated energy barrier for this reaction is about 8 kcal/mol. These predictions await experimental verification. In the case of the neutral ligand TFK0, however, multiple MD simulations on the TFK0/H447I complex reveal that none of the water molecules can be retained in the active site as a "catalytic" water. Furthermore, our alchemical free energy calculation also suggests that the binding of TFK0 to H447I is much weaker than that of TFK+. Taken together, our computational studies confirm that TFK0 is almost inactive in the H447I mutant and also provide detailed mechanistic insights into the experimental observations.
Barriers to Ion Translocation in Cationic and Anionic Receptors from the Cys-Loop FamilyJournal of the American Chemical Society, Vol. 129, Issue 26, pp. 8217-8224 (2007) [PubMed 17552523]
Understanding the mechanisms of gating and ion permeation in biological channels and receptors has been a long-standing challenge in biophysics. Recent advances in structural biology have revealed the architecture of a number of transmembrane channels and allowed detailed, molecular-level insight into these systems. Herein, we have examined the barriers to ion conductance and origins of ion selectivity in models of the cationic human α7 nicotinic acetylcholine receptor (nAChR) and the anionic α1 glycine receptor (GlyR), based on the structure of Torpedo nAChR. Molecular dynamics simulations were used to determine water density profiles along the channel length and established that both receptor pores were fully hydrated. The very low water density in the middle of the nAChR pore indicated the existence of a hydrophobic constriction. By contrast, the pore of GlyR was lined with hydrophilic residues and remained well hydrated throughout. Adaptive biasing force simulations allowed us to reconstruct potentials of mean force (PMFs) for chloride and sodium ions in the two receptors. For the nicotinic receptor we observed barriers to ion translocation associated with rings of hydrophobic residues - Val13' and Leu9' - in the middle of the transmembrane domain. This finding further substantiates the hydrophobic gating hypothesis for nAChR. The PMF revealed no significant hydrophobic barrier for chloride translocation in GlyR. For both receptors non- permeant ions displayed considerable barriers. Thus, the overall electrostatics and the presence of rings of charged residues at the entrance and exit of the channels were sufficient to explain the experimentally observed anion and cation selectivity.
Understanding the mechanisms of gating and ion permeation in biological channels and receptors has been a long-standing challenge in biophysics. Recent advances in structural biology have revealed the architecture of a number of transmembrane channels and allowed detailed, molecular-level insight into these systems. Herein, we have examined the barriers to ion conductance and origins of ion selectivity in models of the cationic human α7 nicotinic acetylcholine receptor (nAChR) and the anionic α1 glycine receptor (GlyR), based on the structure of Torpedo nAChR. Molecular dynamics simulations were used to determine water density profiles along the channel length and established that both receptor pores were fully hydrated. The very low water density in the middle of the nAChR pore indicated the existence of a hydrophobic constriction. By contrast, the pore of GlyR was lined with hydrophilic residues and remained well hydrated throughout. Adaptive biasing force simulations allowed us to reconstruct potentials of mean force (PMFs) for chloride and sodium ions in the two receptors. For the nicotinic receptor we observed barriers to ion translocation associated with rings of hydrophobic residues - Val13' and Leu9' - in the middle of the transmembrane domain. This finding further substantiates the hydrophobic gating hypothesis for nAChR. The PMF revealed no significant hydrophobic barrier for chloride translocation in GlyR. For both receptors non- permeant ions displayed considerable barriers. Thus, the overall electrostatics and the presence of rings of charged residues at the entrance and exit of the channels were sufficient to explain the experimentally observed anion and cation selectivity.
Dealing with bound waters in a site: Do they leave or stay?In "Computational and Structural Approaches to Drug Discovery," R.M. Stroud, Ed., Royal Society of Chemistry, Ch. 5, pp. 95-110 (2007)
Water molecules are ubiquitous to biomolecules. They play a very important role in the function and structure of proteins. They also make up an integral part of protein structures and contribute to their stability. Localized water molecules can be found in the crevices on protein surfaces, in deeper channels or ligand binding sites, and within buried hydrophilic cavities of protein molecules. There is also evidence that water molecules could sometimes bind to cavities that are mostly hydrophobic. Also, water molecules are characteristically present in the interfaces found in protein-ligand, protein-protein, protein-carbohydrate, protein-nucleic acid, and nucleic acid-ligand complexes. These water molecules can now be routinely detected by X-ray crytallographic and nuclear magnetic resonance (NMR) experiments. Interfacial water molecules typically act as bridges between two solute molecules, playing a role in recognition and specificity, and at the same time stabilizing the complex by accepting and donating hydrogen bonds. However, in some cases water molecules have been shown to direct the function of an enzyme protein, for example the catalytic abilities of protein kinases. The highly conserved protein kinase CK2, which utilizes ATP, could also efficiently make use of GTP by localizing a water molecule at the interface between GTP and the binding pocket.
Water molecules are ubiquitous to biomolecules. They play a very important role in the function and structure of proteins. They also make up an integral part of protein structures and contribute to their stability. Localized water molecules can be found in the crevices on protein surfaces, in deeper channels or ligand binding sites, and within buried hydrophilic cavities of protein molecules. There is also evidence that water molecules could sometimes bind to cavities that are mostly hydrophobic. Also, water molecules are characteristically present in the interfaces found in protein-ligand, protein-protein, protein-carbohydrate, protein-nucleic acid, and nucleic acid-ligand complexes. These water molecules can now be routinely detected by X-ray crytallographic and nuclear magnetic resonance (NMR) experiments. Interfacial water molecules typically act as bridges between two solute molecules, playing a role in recognition and specificity, and at the same time stabilizing the complex by accepting and donating hydrogen bonds. However, in some cases water molecules have been shown to direct the function of an enzyme protein, for example the catalytic abilities of protein kinases. The highly conserved protein kinase CK2, which utilizes ATP, could also efficiently make use of GTP by localizing a water molecule at the interface between GTP and the binding pocket.
Remarkable Loop Flexibility in Avian Influenza N1 and its Implications for Antiviral Drug DesignJournal of the American Chemical Society, Vol. 129, Issue 25, pp. 7764-7765 (2007) [PubMed 17539643]
The emergence and continuing global spread of the highly virulent avian influenza H5N1 has raised concerns of a possible human pandemic. Several approved anti-influenza drugs effectively target the neuraminidase (NA), a surface glycoprotein that cleaves terminal sialic acid residues and facilitates the release of viral progeny from infected cells. The first crystal structures of group-1 NAs revealed that although the binding pose of oseltamivir was similar to that seen in previous crystallographic complexes, the 150-loop adopted a distinct conformation, opening a new cavity adjacent to the active site. Here we show that the 150-loop is able to open into significantly wider conformations than seen in the crystal structures, through explicitly solvated MD simulations of the apo and oseltamivir-bound forms of tetrameric N1. We find that motion in the 150-loop is coupled to motion in the neighboring 430-loop, which expands the active site cavity even further. Furthermore, in simulations of the oseltamivir-bound system, the 150-loop approaches the closed conformation, suggesting that the loop switching motion may be more rapid than previously observed.
The emergence and continuing global spread of the highly virulent avian influenza H5N1 has raised concerns of a possible human pandemic. Several approved anti-influenza drugs effectively target the neuraminidase (NA), a surface glycoprotein that cleaves terminal sialic acid residues and facilitates the release of viral progeny from infected cells. The first crystal structures of group-1 NAs revealed that although the binding pose of oseltamivir was similar to that seen in previous crystallographic complexes, the 150-loop adopted a distinct conformation, opening a new cavity adjacent to the active site. Here we show that the 150-loop is able to open into significantly wider conformations than seen in the crystal structures, through explicitly solvated MD simulations of the apo and oseltamivir-bound forms of tetrameric N1. We find that motion in the 150-loop is coupled to motion in the neighboring 430-loop, which expands the active site cavity even further. Furthermore, in simulations of the oseltamivir-bound system, the 150-loop approaches the closed conformation, suggesting that the loop switching motion may be more rapid than previously observed.
“New-version-fast-multipole-method” accelerated electrostatic calculations in biomolecular systemsJournal of Computational Physics, Vol. 226, Issue 2, pp. 1348-1366 (2007) [PubMed 18379638]
In this paper, we present an efficient and accurate numerical algorithm for calculating the electrostatic interactions in biomolecular systems. In our scheme, a boundary integral equation (BIE) approach is applied to discretize the linearized Poisson-Boltzmann (PB) equation. The resulting integral formulas are well conditioned for single molecule cases as well as for systems with more than one macromolecule, and are solved efficiently using Krylov subspace based iterative methods such as generalized minimal residual (GMRES) or biconjugate gradients stabilized (BiCGStab) methods. In each iteration, the convolution type matrix-vector multiplications are accelerated by a new version of the fast multipole method (FMM). The implemented algorithm is asymptotically optimal $O(N)$ both in CPU time and memory usage with optimized prefactors. Our approach enhances the present computational ability to treat electrostatics of large scale systems in protein-protein interactions and nano particle assembly processes. Applications including calculating the electrostatics of the nicotinic acetylcholine receptor (nAChR) and interactions between protein Sso7d and DNA are presented.
In this paper, we present an efficient and accurate numerical algorithm for calculating the electrostatic interactions in biomolecular systems. In our scheme, a boundary integral equation (BIE) approach is applied to discretize the linearized Poisson-Boltzmann (PB) equation. The resulting integral formulas are well conditioned for single molecule cases as well as for systems with more than one macromolecule, and are solved efficiently using Krylov subspace based iterative methods such as generalized minimal residual (GMRES) or biconjugate gradients stabilized (BiCGStab) methods. In each iteration, the convolution type matrix-vector multiplications are accelerated by a new version of the fast multipole method (FMM). The implemented algorithm is asymptotically optimal $O(N)$ both in CPU time and memory usage with optimized prefactors. Our approach enhances the present computational ability to treat electrostatics of large scale systems in protein-protein interactions and nano particle assembly processes. Applications including calculating the electrostatics of the nicotinic acetylcholine receptor (nAChR) and interactions between protein Sso7d and DNA are presented.
Application of the level-set method to the implicit solvation of nonpolar moleculesJournal of Chemical Physics, Vol. 127, article 084503, 10 pages (2007) [PubMed 17764265]
A level-set method is developed for numerically capturing the equilibrium solute-solvent interface that is defined by the recently proposed variational implicit solvent model (http://dx.doi.org/10.1103/PhysRevLett.96.087802" target="_blank" class="ref">Dzubiella, Swanson, and McCammon, Phys. Rev. Lett. 96, 087802 (2006) and http://dx.doi.org/10.1063/1.2171192" target="_blank" class="ref">J. Chem. Phys. 124, 084905 (2006)). In the level-set method, a possible solute-solvent interface is represented by the zero level-set (i.e., the zero level surface) of a level-set function and is eventually evolved into the equilibrium solute-solvent interface. The evolution law is determined by minimization of a solvation free energy functional that couples both the interfacial energy and the van der Waals type solute-solvent interaction energy. The surface evolution is thus an energy minimizing process, and the equilibrium solute-solvent interface is an output of this process. The method is implemented and applied to the solvation of nonpolar molecules such as two xenon atoms, two parallel paraffin plates, helical alkane chains, and a single fullerene C60. The level-set solutions show good agreement for the solvation energy when compared to available molecular dynamics simulations. In particular, the method captures solvent dewetting (nanobubble formation) and quantitatively describes the interaction in the strongly hydrophobic plate system.
A level-set method is developed for numerically capturing the equilibrium solute-solvent interface that is defined by the recently proposed variational implicit solvent model (http://dx.doi.org/10.1103/PhysRevLett.96.087802" target="_blank" class="ref">Dzubiella, Swanson, and McCammon, Phys. Rev. Lett. 96, 087802 (2006) and http://dx.doi.org/10.1063/1.2171192" target="_blank" class="ref">J. Chem. Phys. 124, 084905 (2006)). In the level-set method, a possible solute-solvent interface is represented by the zero level-set (i.e., the zero level surface) of a level-set function and is eventually evolved into the equilibrium solute-solvent interface. The evolution law is determined by minimization of a solvation free energy functional that couples both the interfacial energy and the van der Waals type solute-solvent interaction energy. The surface evolution is thus an energy minimizing process, and the equilibrium solute-solvent interface is an output of this process. The method is implemented and applied to the solvation of nonpolar molecules such as two xenon atoms, two parallel paraffin plates, helical alkane chains, and a single fullerene C60. The level-set solutions show good agreement for the solvation energy when compared to available molecular dynamics simulations. In particular, the method captures solvent dewetting (nanobubble formation) and quantitatively describes the interaction in the strongly hydrophobic plate system.
Nanosecond-Timescale Conformational Dynamics of the Human α7 Nicotinic Acetylcholine ReceptorBiophysical Journal, Vol. 93, Issue 8, pp. 2622-2634 (2007) [PubMed 17573436]
We explore the conformational dynamics of a homology model of the human alpha7 nicotinic acetylcholine receptor (nAChR) using molecular dynamics (MD) simulation and analyses of root-mean-square fluctuations (RMSF), block partitioning of segmental motion and principal component analysis (PCA). The results reveal flexible regions and concerted global motions of the subunits encompassing extracellular and transmembrane domains of the subunits. The most relevant motions comprise a bending, hinged at the beta10-M1 region, accompanied by concerted tilting of the M2 helices that widens the intracellular end of the channel. Despite the nanosecond timescale, the observations suggest that tilting of the M2 helices may initiate opening of the pore. The results also reveal direct coupling between a twisting motion of the extracellular domain and dynamic changes of M2. Covariance analysis of inter-residue motions shows that this coupling arises through a network of residues within the Cys- and M2-M3 loops where Phe135 is stabilized within a hydrophobic pocket formed by Leu270 and Ile271. The resulting concerted motion causes a downward shift of the M2 helices that disrupts a hydrophobic girdle formed by 9' and 13' residues.
We explore the conformational dynamics of a homology model of the human alpha7 nicotinic acetylcholine receptor (nAChR) using molecular dynamics (MD) simulation and analyses of root-mean-square fluctuations (RMSF), block partitioning of segmental motion and principal component analysis (PCA). The results reveal flexible regions and concerted global motions of the subunits encompassing extracellular and transmembrane domains of the subunits. The most relevant motions comprise a bending, hinged at the beta10-M1 region, accompanied by concerted tilting of the M2 helices that widens the intracellular end of the channel. Despite the nanosecond timescale, the observations suggest that tilting of the M2 helices may initiate opening of the pore. The results also reveal direct coupling between a twisting motion of the extracellular domain and dynamic changes of M2. Covariance analysis of inter-residue motions shows that this coupling arises through a network of residues within the Cys- and M2-M3 loops where Phe135 is stabilized within a hydrophobic pocket formed by Leu270 and Ile271. The resulting concerted motion causes a downward shift of the M2 helices that disrupts a hydrophobic girdle formed by 9' and 13' residues.
Functional and Structural Insights Revealed by Molecular Dynamics Simulations of an Essential RNA Editing Ligase in Trypanosoma bruceiPLoS Neglected Tropical Diseases, Vol. 1, No. 2, article e68, 10 pages (2007) [PubMed 18060084]
RNA editing ligase 1 (TbREL1) is required for the survival of both the insect and bloodstream forms of Trypanosoma brucei, the parasite responsible for the devastating tropical disease, African sleeping sickness. The type of RNA editing that TbREL1 is involved in is unique to the trypanosomes, and no close human homolog is known to exist. In addition, the high-resolution crystal structure revealed several unique features of the active site, making this enzyme a promising target for structure-based drug design. In this work, two 20 ns atomistic molecular dynamics (MD) simulations are employed to investigate the dynamics of TbREL1, both with and without the ATP substrate present. The flexibility of the active site, dynamics of conserved residues and crystallized water molecules, and the interactions between TbREL1 and the ATP substrate are investigated and discussed in the context of TbREL1's function. Differences in local and global motion upon ATP binding suggest that two peripheral loops, unique to the trypanosomes, may be involved in interdomain signaling events. Notably, a significant structural rearrangement of the enzyme's active site occurs during the apo simulations, opening an additional cavity adjacent to the ATP binding site that could be exploited in the development of effective inhibitors directed against this protozoan parasite. Finally, ensemble averaged electrostatics calculations over the MD simulations reveal a novel putative RNA binding site, a discovery that has previously eluded scientists. Ultimately, we use the insights gained through the MD simulations to make several predictions and recommendations, which we anticipate will help direct future experimental studies and structure-based drug discovery efforts against this vital enzyme.
RNA editing ligase 1 (TbREL1) is required for the survival of both the insect and bloodstream forms of Trypanosoma brucei, the parasite responsible for the devastating tropical disease, African sleeping sickness. The type of RNA editing that TbREL1 is involved in is unique to the trypanosomes, and no close human homolog is known to exist. In addition, the high-resolution crystal structure revealed several unique features of the active site, making this enzyme a promising target for structure-based drug design. In this work, two 20 ns atomistic molecular dynamics (MD) simulations are employed to investigate the dynamics of TbREL1, both with and without the ATP substrate present. The flexibility of the active site, dynamics of conserved residues and crystallized water molecules, and the interactions between TbREL1 and the ATP substrate are investigated and discussed in the context of TbREL1's function. Differences in local and global motion upon ATP binding suggest that two peripheral loops, unique to the trypanosomes, may be involved in interdomain signaling events. Notably, a significant structural rearrangement of the enzyme's active site occurs during the apo simulations, opening an additional cavity adjacent to the ATP binding site that could be exploited in the development of effective inhibitors directed against this protozoan parasite. Finally, ensemble averaged electrostatics calculations over the MD simulations reveal a novel putative RNA binding site, a discovery that has previously eluded scientists. Ultimately, we use the insights gained through the MD simulations to make several predictions and recommendations, which we anticipate will help direct future experimental studies and structure-based drug discovery efforts against this vital enzyme.
Free Energy Profile of H-ras Membrane Anchor upon Membrane InsertionAngewandte Chemie International Edition, Vol. 46, Issue 43, pp. 8234-8237 (2007) [PubMed 17886310]
Ras GTPases mediatr signaling pathways in cell proliferation, development, and apoptosis. They undergo isoprenylation at a C-terminal CaaX signal (a usually represents aliphatic and X any amino acid) followed by proteolysis of aaX and carboxymethylation. In the case of H-ras, a subsequent dual palmitoylation of cysteines adjacent to the site of farnesylation produces a mature anchor for plasma membrane targeting. http://dx.doi.org/10.1038/nchembio834" target="_blank" class="ref">M.D. Resh, Nat. Chem. Biol. 2, 584 (2006); http://www.pubmedcentral.nih.gov/articlerender.fcgi?pubmedid=2663468" target="_blank" class="ref">L. Gutierrez, A.I. Magee, C.J. Marshall, J.F. Hancock, EMBO J. 8, 1093 (1989) Atomistic information, such as the structure of membrane-bound ras and the free energy of complex formation, are vital in research efforts geared towards designing ras-isoform-selective anticancer agents. The most common experimental techniques are not yet able to provide such information. Here we present computational results on the free energy profile for the transfer of the H-ras membrane anchor from water to a bilayer of 1,2-dimyristoyl-sn-glycero-3-phosphocholine (1,2-DMPC) lipids. We find that there is no significant barrier for insertion, and that once a few carbon atoms of the ras lipid chains cross the membrane-water interface, the free energy displays a steeply downhill profile. Insertion into the hydrocarbon core of the ras lipids and the interfacial localization of the backbone together produce a gain in free energy of up to 30 kcal/mol. Additionally, using the recently reported computationally derived structures of full-length H-ras in a DMPC bilayer, http://dx.doi.org/10.1021/jm061053f" target="_blank" class="ref">A.A. Gorfe, M. Hanzal-Bayer, D. Abankwa, J.F. Hancock, J.A. McCammon, J. Med. Chem. 50, 674 (2007) we explain how a small difference in free energy would enable modulation of H-ras membrane binding by the linker and the catalytic domain.
Ras GTPases mediatr signaling pathways in cell proliferation, development, and apoptosis. They undergo isoprenylation at a C-terminal CaaX signal (a usually represents aliphatic and X any amino acid) followed by proteolysis of aaX and carboxymethylation. In the case of H-ras, a subsequent dual palmitoylation of cysteines adjacent to the site of farnesylation produces a mature anchor for plasma membrane targeting. http://dx.doi.org/10.1038/nchembio834" target="_blank" class="ref">M.D. Resh, Nat. Chem. Biol. 2, 584 (2006); http://www.pubmedcentral.nih.gov/articlerender.fcgi?pubmedid=2663468" target="_blank" class="ref">L. Gutierrez, A.I. Magee, C.J. Marshall, J.F. Hancock, EMBO J. 8, 1093 (1989) Atomistic information, such as the structure of membrane-bound ras and the free energy of complex formation, are vital in research efforts geared towards designing ras-isoform-selective anticancer agents. The most common experimental techniques are not yet able to provide such information. Here we present computational results on the free energy profile for the transfer of the H-ras membrane anchor from water to a bilayer of 1,2-dimyristoyl-sn-glycero-3-phosphocholine (1,2-DMPC) lipids. We find that there is no significant barrier for insertion, and that once a few carbon atoms of the ras lipid chains cross the membrane-water interface, the free energy displays a steeply downhill profile. Insertion into the hydrocarbon core of the ras lipids and the interfacial localization of the backbone together produce a gain in free energy of up to 30 kcal/mol. Additionally, using the recently reported computationally derived structures of full-length H-ras in a DMPC bilayer, http://dx.doi.org/10.1021/jm061053f" target="_blank" class="ref">A.A. Gorfe, M. Hanzal-Bayer, D. Abankwa, J.F. Hancock, J.A. McCammon, J. Med. Chem. 50, 674 (2007) we explain how a small difference in free energy would enable modulation of H-ras membrane binding by the linker and the catalytic domain.
Generalized gradient-augmented harmonic Fourier beads method with multiple atomic and/or center-of-mass positional restraintsJournal of Chemical Physics, Vol. 127, article 124901, 12 pages (2007) [PubMed 17902931]
We describe a generalization of the gradient-augmented Harmonic Fourier Beads method for finding minimum free energy transition path ensembles and similarly minimum potential energy paths to allow positional restraints on the centers of mass of selected atoms. The generalized method (ggaHFB) further extends the scope of the HFB methodology to studying molecule transport across various mobile phases such as lipid membranes. Furthermore, the new implementation improves the applicability of the HFB method to studies of ligand binding, protein folding and enzyme catalysis as well as modeling equilibrium pulling experiments. Like its predecessor, the ggaHFB method provides accurate energy profiles along the specified paths and in certain simple cases avoids the need for path optimization. The utility of the ggaHFB method is demonstrated with an application to the water permeation through a single-wall (5,5) carbon nanotube with a diameter of 6.78 A and length of 16.0 A. We provide a simple rationale as to why water enters the hydrophobic nanotube and why it does so in pulses and in wire assembly.
We describe a generalization of the gradient-augmented Harmonic Fourier Beads method for finding minimum free energy transition path ensembles and similarly minimum potential energy paths to allow positional restraints on the centers of mass of selected atoms. The generalized method (ggaHFB) further extends the scope of the HFB methodology to studying molecule transport across various mobile phases such as lipid membranes. Furthermore, the new implementation improves the applicability of the HFB method to studies of ligand binding, protein folding and enzyme catalysis as well as modeling equilibrium pulling experiments. Like its predecessor, the ggaHFB method provides accurate energy profiles along the specified paths and in certain simple cases avoids the need for path optimization. The utility of the ggaHFB method is demonstrated with an application to the water permeation through a single-wall (5,5) carbon nanotube with a diameter of 6.78 A and length of 16.0 A. We provide a simple rationale as to why water enters the hydrophobic nanotube and why it does so in pulses and in wire assembly.
Solvent reaction field potential inside an uncharged globular protein: A bridge between implicit and explicit solvent models?Journal of Chemical Physics, Vol. 127, article 155101, 12 pages (2007) [PubMed 17949217]
The solvent reaction field potential of an uncharged protein immersed in Simple Point Charge / Extended (SPC/E) explicit solvent was computed over a series of molecular dynamics trajectories, in total 1560 ns of simulation time. A finite, positive potential of 13 to 24 kT/e (where T = 300K), dependent on the geometry of the solvent-accessible surface, was observed inside the biomolecule. The primary contribution to this potential arose from a layer of positive charge density 1.0 Angstroms from the solute surface, on average 0.008 e per cubic Angstrom, which we found to be the product of a highly ordered first solvation shell. Significant second solvation shell effects, including additional layers of charge density and a slight decrease in the short-range solvent-solvent interaction strength, were also observed. The impact of these findings on implicit solvent models was assessed by running similar explicit-solvent simulations on the fully charged protein system. When the energy due to the solvent reaction field in the uncharged system is accounted for, correlation between per-atom electrostatic energies for the explicit solvent model and a simple implicit (Poisson) calculation is 0.97, and correlation between per-atom energies for the explicit solvent model and a previously published, optimized Poisson model is 0.99.
The solvent reaction field potential of an uncharged protein immersed in Simple Point Charge / Extended (SPC/E) explicit solvent was computed over a series of molecular dynamics trajectories, in total 1560 ns of simulation time. A finite, positive potential of 13 to 24 kT/e (where T = 300K), dependent on the geometry of the solvent-accessible surface, was observed inside the biomolecule. The primary contribution to this potential arose from a layer of positive charge density 1.0 Angstroms from the solute surface, on average 0.008 e per cubic Angstrom, which we found to be the product of a highly ordered first solvation shell. Significant second solvation shell effects, including additional layers of charge density and a slight decrease in the short-range solvent-solvent interaction strength, were also observed. The impact of these findings on implicit solvent models was assessed by running similar explicit-solvent simulations on the fully charged protein system. When the energy due to the solvent reaction field in the uncharged system is accounted for, correlation between per-atom electrostatic energies for the explicit solvent model and a simple implicit (Poisson) calculation is 0.97, and correlation between per-atom energies for the explicit solvent model and a previously published, optimized Poisson model is 0.99.
H-ras Protein in a Bilayer: Interaction and Structure PerturbationJournal of the American Chemical Society, Vol. 129, No. 40, pp. 12280-12286 (2007) [PubMed 17880077]
Ras GTPases become functionally active when anchored to membranes by inserting their lipid modified side chains. Their role in cell division, development and cancer has made them targets of extensive research efforts, yet the mechanism of membrane insertion and the structure of the resulting complex remain elusive. Recently, the structure of the full-length H-ras protein in a DMPC bilayer has been computationally characterized. Here, the atomic interactions between the H-ras membrane anchor and the DMPC bilayer are investigated in detail. We find that the palmitoylated cysteines and Met182 have dual contributions to membrane affinity: hydrogen bonding by their amides and vdW interaction by their hydrophobic side chains. The polar side chains help maintain the orientation of the anchor. Although the overall structure of the bilayer is similar to that of a pure DMPC, there are localized perturbations. These perturbations depend on the insertion depth and backbone localization of the anchor, which in turn is modulated by the catalytic domain and the linker. The pattern of anchor amide-DMPC phosphate/carbonyl hydrogen bonds and the flexibility of Palm184 are important in discriminating between different modes of ras-DMPC interactions. The results provide structural arguments in support of the proposed participation of ras in the organization of membrane nanoclusters.
Ras GTPases become functionally active when anchored to membranes by inserting their lipid modified side chains. Their role in cell division, development and cancer has made them targets of extensive research efforts, yet the mechanism of membrane insertion and the structure of the resulting complex remain elusive. Recently, the structure of the full-length H-ras protein in a DMPC bilayer has been computationally characterized. Here, the atomic interactions between the H-ras membrane anchor and the DMPC bilayer are investigated in detail. We find that the palmitoylated cysteines and Met182 have dual contributions to membrane affinity: hydrogen bonding by their amides and vdW interaction by their hydrophobic side chains. The polar side chains help maintain the orientation of the anchor. Although the overall structure of the bilayer is similar to that of a pure DMPC, there are localized perturbations. These perturbations depend on the insertion depth and backbone localization of the anchor, which in turn is modulated by the catalytic domain and the linker. The pattern of anchor amide-DMPC phosphate/carbonyl hydrogen bonds and the flexibility of Palm184 are important in discriminating between different modes of ras-DMPC interactions. The results provide structural arguments in support of the proposed participation of ras in the organization of membrane nanoclusters.
A proposed signaling motif for nuclear import in mRNA processing via the formation of Arginine ClawProceedings of the National Academy of Sciences of the USA, Vol. 104, No. 38, pp. 14947-14951 (2007) [PubMed 17823247]
Phosphorylation of proteins by kinases is the most commonly studied class of posttranslational modification, yet its structural consequences are not well understood. The human SR (serine-arginine) protein ASF/SF2 relies on the processive phosphorylation of the serine residues of eight consecutive arginine-serine (RS) dipeptide repeats at the C terminus by SRPK1 before it can be transported into the nucleus. This SR protein plays critical roles in spliceosome assembly, pre-mRNA splicing, and mRNA export, and the phosphorylation process of the RS repeats has been extensively studied experimentally. However, knowledge of the conformational changes associated with the phosphorylation of this simple sequence and how it triggers the importation of the SR protein is lacking. Here, we have carried out extensive molecular dynamics simulations to show that phosphorylation of the eight RS repeats significantly alters the peptide's conformation and leads to the formation of very stable structures that are likely to be involved in the recognition, binding, and transport of the SR protein. Specifically, we found an unusual symmetry-broken phase of conformations of the repetitive and quasi-symmetric phosphorylated peptide sequence. One of the main characteristics of these conformations is the exposed phosphate groups on the periphery, which possibly could serve as the recognition platform for the transport protein transportin-SR2.
Phosphorylation of proteins by kinases is the most commonly studied class of posttranslational modification, yet its structural consequences are not well understood. The human SR (serine-arginine) protein ASF/SF2 relies on the processive phosphorylation of the serine residues of eight consecutive arginine-serine (RS) dipeptide repeats at the C terminus by SRPK1 before it can be transported into the nucleus. This SR protein plays critical roles in spliceosome assembly, pre-mRNA splicing, and mRNA export, and the phosphorylation process of the RS repeats has been extensively studied experimentally. However, knowledge of the conformational changes associated with the phosphorylation of this simple sequence and how it triggers the importation of the SR protein is lacking. Here, we have carried out extensive molecular dynamics simulations to show that phosphorylation of the eight RS repeats significantly alters the peptide's conformation and leads to the formation of very stable structures that are likely to be involved in the recognition, binding, and transport of the SR protein. Specifically, we found an unusual symmetry-broken phase of conformations of the repetitive and quasi-symmetric phosphorylated peptide sequence. One of the main characteristics of these conformations is the exposed phosphate groups on the periphery, which possibly could serve as the recognition platform for the transport protein transportin-SR2.
Electrodiffusion: A continuum modeling framework for biomolecular systems with realistic spatiotemporal resolutionJournal of Chemical Physics, Vol. 127, article 135102, 17 pages (2007) [PubMed 17919055]
A computational framework is presented for the continuum modeling of cellular biomolecular diffusion influenced by electrostatic driving forces. This framework is developed from a combination of state-of-the-art numerical methods, geometric meshing and computer visualization tools. In particular, a hybrid of (adaptive) finite element and boundary element methods is adopted to solve the Smoluchowski equation (SE), the Poisson equation (PE), and the Poisson-Nernst-Planck equation (PNPE) in order to describe electrodiffusion processes. The finite element method is used because of its flexibility in modeling irregular geometries and complex boundary conditions. The boundary element method is used due to the convenience of treating the singularities in the source charge distribution and its accurate solution to electrostatic problems on molecular boundaries. Nonsteady-state diffusion can be studied using this framework, with the electric field computed using the densities of charged small molecules and mobile ions in the solvent. A solution for mesh generation for biomolecular systems is supplied, which is an essential component for the finite element and boundary element computations. The uncoupled Smoluchowski equation and Poisson-Boltzmann equation (PBE) are considered as special cases of the PNPE in numerical algorithm, and therefore can be solved in this framework as well. Two types of computations are reported in the results: stationary PNPE and time-dependent SE or Nernst-Planck equations (PN) solutions. A biological application of the first type is the ionic density distribution around a fragment of DNA determined by the equilibrium PNPE. The stationary PNPE with non-zero flux is also discussed for a simple model system. The second is a time-dependent diffusion process: the consumption of the neurotransmitter acetylcholine (ACh) by acetylcholinesterase (AChE), determined by the SE and a single uncoupled solution of the Poisson-Boltzmann equation. The electrostatic effects, counterion compensation, spatiotemporal distribution, and diffusion-controlled reaction kinetics are analyzed and different methods are compared.
A computational framework is presented for the continuum modeling of cellular biomolecular diffusion influenced by electrostatic driving forces. This framework is developed from a combination of state-of-the-art numerical methods, geometric meshing and computer visualization tools. In particular, a hybrid of (adaptive) finite element and boundary element methods is adopted to solve the Smoluchowski equation (SE), the Poisson equation (PE), and the Poisson-Nernst-Planck equation (PNPE) in order to describe electrodiffusion processes. The finite element method is used because of its flexibility in modeling irregular geometries and complex boundary conditions. The boundary element method is used due to the convenience of treating the singularities in the source charge distribution and its accurate solution to electrostatic problems on molecular boundaries. Nonsteady-state diffusion can be studied using this framework, with the electric field computed using the densities of charged small molecules and mobile ions in the solvent. A solution for mesh generation for biomolecular systems is supplied, which is an essential component for the finite element and boundary element computations. The uncoupled Smoluchowski equation and Poisson-Boltzmann equation (PBE) are considered as special cases of the PNPE in numerical algorithm, and therefore can be solved in this framework as well. Two types of computations are reported in the results: stationary PNPE and time-dependent SE or Nernst-Planck equations (PN) solutions. A biological application of the first type is the ionic density distribution around a fragment of DNA determined by the equilibrium PNPE. The stationary PNPE with non-zero flux is also discussed for a simple model system. The second is a time-dependent diffusion process: the consumption of the neurotransmitter acetylcholine (ACh) by acetylcholinesterase (AChE), determined by the SE and a single uncoupled solution of the Poisson-Boltzmann equation. The electrostatic effects, counterion compensation, spatiotemporal distribution, and diffusion-controlled reaction kinetics are analyzed and different methods are compared.
Dynamics, Hydration, and Motional Averaging of a Loop-Gated Artificial Protein Cavity: The W191G Mutant of Cytochrome c Peroxidase in Water as Revealed by Molecular Dynamics SimulationsBiochemistry, Vol. 46, Issue 37, pp. 10629-10642 (2007) [PubMed 17718514]
Five molecular dynamics simulations of the W191G cavity mutant of cytochrome c peroxidase in explicit water reveal distinct dynamic and hydration behavior depending on the closed or open state of the flexible loop gating the cavity, the binding of (K+ or small molecule) cations, and the system temperature. The conformational spaces sampled by the loop region and by the cavity significantly reduce upon binding. The largest ordering factor on water dynamics is the presence of the K+ ion occupying the gated cavity. Considerable water exchange occurs for the open-gate cavity when no ligand or cation is bound. In all cases, good correspondence is found between the calculated (ensemble-averaged) location of water molecules and the water sites determined by X-ray crystallography experiments. However, our simulations suggest that these sites do not necessarily correspond to the presence of bound water molecules. In fact, individual water molecules may repeatedly exchange within the cavity volume yet occupy on average these water sites. Four major conclusions emerge. First, it seems misleading to interpret the conformation of protein loop regions in terms of single dominant structures. Second, our simulations support the general picture of Pro 190 cis-trans isomerization as a determinant of the loop-opening mechanism. Third, receptor flexibility is fundamental for ligand binding and molecular recognition, and our results suggest its importance for the docking of small compounds to the artificial cavity. Fourth, after validation against the available experimental data, molecular dynamics simulations can be used to characterize the dynamics and exchange of water molecules and ions, providing atomic level and time- dependent information otherwise inaccessible to experiments.
Five molecular dynamics simulations of the W191G cavity mutant of cytochrome c peroxidase in explicit water reveal distinct dynamic and hydration behavior depending on the closed or open state of the flexible loop gating the cavity, the binding of (K+ or small molecule) cations, and the system temperature. The conformational spaces sampled by the loop region and by the cavity significantly reduce upon binding. The largest ordering factor on water dynamics is the presence of the K+ ion occupying the gated cavity. Considerable water exchange occurs for the open-gate cavity when no ligand or cation is bound. In all cases, good correspondence is found between the calculated (ensemble-averaged) location of water molecules and the water sites determined by X-ray crystallography experiments. However, our simulations suggest that these sites do not necessarily correspond to the presence of bound water molecules. In fact, individual water molecules may repeatedly exchange within the cavity volume yet occupy on average these water sites. Four major conclusions emerge. First, it seems misleading to interpret the conformation of protein loop regions in terms of single dominant structures. Second, our simulations support the general picture of Pro 190 cis-trans isomerization as a determinant of the loop-opening mechanism. Third, receptor flexibility is fundamental for ligand binding and molecular recognition, and our results suggest its importance for the docking of small compounds to the artificial cavity. Fourth, after validation against the available experimental data, molecular dynamics simulations can be used to characterize the dynamics and exchange of water molecules and ions, providing atomic level and time- dependent information otherwise inaccessible to experiments.
Sampling of slow diffusive conformational transitions with accelerated molecular dynamicsJournal of Chemical Physics, Vol. 127, article 155102, 9 pages (2007) [PubMed 17949218]
Slow diffusive conformational transitions play key functional roles in biomolecular systems. Our ability to sample these motions with molecular dynamics simulation in explicit solvent is limited by the slow diffusion of the solvent molecules around the biomolecules. Previously, we proposed an accelerated molecular dynamics method that has been shown to efficiently sample the torsional degrees of freedom of biomolecules beyond the millisecond timescale. However, in our previous approach, large-amplitude displacements of biomolecules are still slowed by the diffusion of the solvent. Here we present a unified approach of efficiently sampling both the torsional degrees of freedom and the diffusive motions concurrently. We show that this approach samples the configuration space more efficiently than normal molecular dynamics and that ensemble averages converge faster to the correct values.
Slow diffusive conformational transitions play key functional roles in biomolecular systems. Our ability to sample these motions with molecular dynamics simulation in explicit solvent is limited by the slow diffusion of the solvent molecules around the biomolecules. Previously, we proposed an accelerated molecular dynamics method that has been shown to efficiently sample the torsional degrees of freedom of biomolecules beyond the millisecond timescale. However, in our previous approach, large-amplitude displacements of biomolecules are still slowed by the diffusion of the solvent. Here we present a unified approach of efficiently sampling both the torsional degrees of freedom and the diffusive motions concurrently. We show that this approach samples the configuration space more efficiently than normal molecular dynamics and that ensemble averages converge faster to the correct values.
Estimating Kinetic Rates from Accelerated Molecular Dynamics Simulations: Alanine Dipeptide in Explicit Solvent as a Case StudyJournal of Chemical Physics, Vol. 127, article 175105, 8 pages (2007) [PubMed 17994855]
Molecular dynamics (MD) simulation is the standard computational technique used to obtain information on the time evolution of the conformations of proteins and many other molecular systems. However, for most biological systems of interest, the time scale for slow conformational transitions is still inaccessible to standard MD simulations. Several sampling methods have been proposed to address this issue, including the accelerated molecular dynamics method. In this work, we study the extent of sampling of the phi/psi space of alanine dipeptide in explicit water using accelerated molecular dynamics and present a framework to recover the correct kinetic rate constant for the helix to beta-strand transition. We show that the accelerated MD can drastically enhance the sampling of the phi/psi conformational phase space when compared to normal MD. In addition, the free energy density plots of the phi/psi space show that all minima regions are accurately sampled and the canonical distribution is recovered. Moreover, the kinetic rate constant for the helix to beta-strand transition is accurately estimated from these simulations by relating the diffusion coefficient to the local energetic roughness of the energy landscape. Surprisingly, even for such a low barrier transition, it is difficult to obtain enough transitions to accurately estimate the rate constant when one uses normal MD.
Molecular dynamics (MD) simulation is the standard computational technique used to obtain information on the time evolution of the conformations of proteins and many other molecular systems. However, for most biological systems of interest, the time scale for slow conformational transitions is still inaccessible to standard MD simulations. Several sampling methods have been proposed to address this issue, including the accelerated molecular dynamics method. In this work, we study the extent of sampling of the phi/psi space of alanine dipeptide in explicit water using accelerated molecular dynamics and present a framework to recover the correct kinetic rate constant for the helix to beta-strand transition. We show that the accelerated MD can drastically enhance the sampling of the phi/psi conformational phase space when compared to normal MD. In addition, the free energy density plots of the phi/psi space show that all minima regions are accurately sampled and the canonical distribution is recovered. Moreover, the kinetic rate constant for the helix to beta-strand transition is accurately estimated from these simulations by relating the diffusion coefficient to the local energetic roughness of the energy landscape. Surprisingly, even for such a low barrier transition, it is difficult to obtain enough transitions to accurately estimate the rate constant when one uses normal MD.
Accelerated entropy estimates with accelerated dynamicsJournal of Chemical Physics, Vol. 127, article 154105, 5 pages (2007) [PubMed 17949130]
Evaluating the entropy of a classical ensemble, requires integrating over the probability density, p(r), as a function of configurational space position, r. In practice, the actual probability is rarely known and the entropy is calculated using a density estimate based on statistical sampling. In complex systems such as biological macromolecules, sampling all significant configurations is notoriously difficult. Although technological progress has enabled molecular dynamics simulations, which integrate Newton's equations of motion to observe time-dependent system behaviors, to reach the nano- to microsecond regime, some biologically relevant conformational changes are known to occur on the order of seconds, hours, or even days. Even the longest trajectories can be kinetically trapped in local energy minima and fall far short of ergodicity. Monte Carlo algorithms, in which random moves are proposed and stochastically accepted according to their energetic properties, encounter a similar problem: it is difficult to design large moves in phase space with reasonable acceptance probabilities.
Evaluating the entropy of a classical ensemble, requires integrating over the probability density, p(r), as a function of configurational space position, r. In practice, the actual probability is rarely known and the entropy is calculated using a density estimate based on statistical sampling. In complex systems such as biological macromolecules, sampling all significant configurations is notoriously difficult. Although technological progress has enabled molecular dynamics simulations, which integrate Newton's equations of motion to observe time-dependent system behaviors, to reach the nano- to microsecond regime, some biologically relevant conformational changes are known to occur on the order of seconds, hours, or even days. Even the longest trajectories can be kinetically trapped in local energy minima and fall far short of ergodicity. Monte Carlo algorithms, in which random moves are proposed and stochastically accepted according to their energetic properties, encounter a similar problem: it is difficult to design large moves in phase space with reasonable acceptance probabilities.
