Chaos in protein dynamics

Stochastic Lag Time Parameterization for Markov State Models of Protein Dynamics

Markov state models (MSMs) play a key role in studying protein conformational dynamics. A sliding count window with a fixed lag time is widely used to sample sub-trajectories for transition counting and MSM construction. However, sub-trajectories sampled with a fixed lag time may not perform well under different selections of lag time, which requires strong prior practice and leads to less robust estimation. To alleviate it, we propose a novel stochastic method from a Poisson process to generate perturbative lag time for sub-trajectory sampling and utilize it to construct a Markov chain. Comprehensive evaluations on the double-well system, WW domain, BPTI, and RBD-ACE2 complex of SARS-CoV-2 reveal that our algorithm significantly increases the robustness and power of a constructed MSM without disturbing the Markovian properties. Furthermore, the superiority of our algorithm is amplified for slow dynamic modes in complex biological processes.

Physicochemical Insights on Terahertz Wave Diminished Side Effects of Drugs from Slow Dissociation

Dopamine D2 receptors (D2Rs) are one of the most intensely investigated and well-established drug targets for neuropsychiatric disorders. Selective D2R antagonists have been developed as efficacious antipsychotic drugs. Nevertheless, the potent drugs with necessarily high affinity are prone to slow dissociation, which invokes a plethora of severe side effects such as extrapyramidal symptoms, substantial weight gain, associated diabetes, etc. This has become a major barrier in treating psychiatric patients. In this work, we propose a physical method, terahertz wave modulation, to promote the dissociation of high-affinity antipsychotics and thus diminish the side effects. We have proven that a 4.0 THz wave could reduce the affinity by 71% between the D2R and a risperidone ligand and meanwhile expand the exit via conformation modulation, which promises an accelerated dissociation of risperidone. In addition, it is estimated that the enhancement of the dissociation rate due to lowered binding by terahertz irritation could constitute up to 8 orders of magnitude, which is fairly impressive and resembles the enzyme's catalysis. Also, acceleration of the dissociation rate could be adjusted by the irritation strength. This work elaborates the terahertz wave-modulated noncovalent interactions critical in cell signaling pathways. Most importantly, it demonstrates the feasibility of terahertz technologies intervening in receptor-ligand complex regulated diseases such as neurodegenerative disorders, metabolic diseases, etc.

Deterministic chaos in the self-assembly of β sheet nanotubes from an amphipathic oligopeptide

The self-assembly of designed peptides into filaments and other higher-order structures has been the focus of intense interest because of the potential for creating new biomaterials and biomedical devices. These peptide assemblies have also been used as models for understanding biological processes, such as the pathological formation of amyloid. We investigate the assembly of an octapeptide sequence, Ac-FKFEFKFE-NH₂, motivated by prior studies that demonstrated that this amphipathic β strand peptide self-assembled into fibrils and biocompatible hydrogels. Using high-resolution cryoelectron microscopy (cryo-EM), we are able to determine the atomic structure for two different coexisting forms of the fibrils, containing four and five β sandwich protofilaments, respectively. Surprisingly, the inner walls in both forms are parallel β sheets, while the outer walls are antiparallel β sheets. Our results demonstrate the chaotic nature of peptide self-assembly and illustrate the importance of cryo-EM structural analysis to understand the complex phase behavior of these materials at near-atomic resolution.

Molecular Insights and Functional Consequences of the Interaction of Heme with Activated Protein C

Aims: In hemolysis, which is accompanied by increased levels of labile redox-active heme and is often associated with hemostatic abnormalities, a decreased activity of activated protein C (APC) is routinely detected. APC is a versatile enzyme that exerts its anticoagulant function through inactivation of clotting factors Va and VIIIa. APC has not been demonstrated to be affected by heme as described for other clotting factors and, thus, is a subject of investigation. Results: We report the interaction of heme with APC and its impact on the protein function by employing spectroscopic and physiologically relevant methods. Binding of heme to APC results in inhibition of its amidolytic and anticoagulant activity, increase of the peroxidase-like activity of heme, and protection of human umbilical vein endothelial cells from heme-induced hyperpermeability. To define the sites that are responsible for heme binding, we mapped the surface of APC for potential heme-binding motifs. T²⁸⁵GWGYHSSR²⁹³ and W³⁸⁷IHGHIRDK³⁹⁵, both located on the basic exosite, turned out as potential heme-binding sites. Molecular docking employing a homology model of full-length APC indicated Tyr²⁸⁹ and His³⁹¹ as the Fe(III)-coordinating amino acids. Innovation: The results strongly suggest that hemolysis-derived heme may directly influence the protein C pathway through binding to APC, conceivably explaining the decreased activity of APC under hemolytic conditions. Further, these results extend our understanding of heme as a multifaceted effector molecule within coagulation and may allow for an improved understanding of disease development in hemostasis under hemolytic conditions. Conclusion: Our study identifies APC as a heme-binding protein and provides insights into the functional consequences.

Kinetic coherence underlies the dynamics of disordered proteins

The dynamics of two proteins of similar size, the globular lysozyme and the intrinsically disordered Huntingtin interacting protein, has been simulated in three states resembling a globule, a pre-molten globule, and a molten globule. A coherence time τ has been defined, measuring the delay in the display of a stochastic behaviour after a perturbation of the system. This time has been computed for two sets of collective variables: the projection of the phase point onto the positions and momenta subspaces (τ_r and τ_p), and the principal components (PCs) of positions q and momenta π produced by a covariance analysis in these subspaces (τ_q and τ_π). In all states τ_p ≈ 3.5τ_r, and τ_π ≈ 3.5τ_q. The coherence times of individual PCs, τ^(l)_q and τ^(l)_π, have also been computed, and τ^(l)_π > τ^(l)_q in all states. The prevalence of τ_p over τ_r, or of τ_π over τ_q, drives the dynamics of the protein over a time range of ≈1–2 ps; moreover, a hidden synchronism appears to raise the momenta subspace's coherence above that of its individual PCs. In the transition of lysozyme to the molten globule the τ^(l)_q decrease but, unexpectedly, the τ^(l)_π increase; after this transition τ_p ≈ 5τ_r and τ_π ≈ 5τ_q. A gain of kinetic coherence accompanies thus the loss of structural coherence caused by the denaturation of the protein in the transition from globule to molten globule. The increase of the τ^(l)_π does not take place in the analogous transition of the Huntingtin protein. These results are compared with those of a similar analysis performed on three pseudo-proteins designed by scrambling the primary sequence of the Huntingtin interacting protein, and on two oligopeptides. The hidden synchronism appears to be a generic property of these polypeptides. The τ^(l)_π spectrum is similar in denaturated and in intrinsically disordered biomolecules; but the gain of kinetic coherence as a result of denaturation seems to be a specific property of the biologically functional lysozyme.

In the phase space of a globular or intrinsically disordered protein, the momenta's dynamics is less chaotic than the coordinates' dynamics. When a protein is denaturated, a gain in kinetic coherence accompanies the loss of structural coherence.

Convergence versus Divergence Behaviors of Asynchronous Iterations, and Their Applications in Concrete Situations

Asynchronous iterations have long been used in distributed computing algorithms to produce calculation methods that are potentially faster than a serial or parallel approach, but whose convergence is more difficult to demonstrate. Conversely, over the past decade, the study of the complex dynamics of asynchronous iterations has been initiated and deepened, as well as their use in computer security and bioinformatics. The first work of these studies focused on chaotic discrete dynamical systems, and links were established between these dynamics on the one hand, and between random or complex behaviours in the sense of the theory of the same name. Computer security applications have focused on pseudo-random number generation, hash functions, hidden information, and various security aspects of wireless sensor networks. At the bioinformatics level, this study of complex systems has allowed an original approach to understanding the evolution of genomes and protein folding. These various contributions are detailed in this review article, which is an extension of the paper “An update on the topological properties of asynchronous iterations” presented during the Sixth International Conference on Parallel, Distributed, GPU and Cloud Computing (Pareng 2019).

Minimalist Protocell Design: A Molecular System Based Solely on Proteins that Form Dynamic Vesicular Membranes Embedding Enzymatic Functions

Life in its molecular context is characterized by the challenge of orchestrating structure, energy and information processes through compartmentalization and chemical transformations amenable to mimicry of protocell models. Here we present an alternative protocell model incorporating dynamic membranes based on amphiphilic elastin-like proteins (ELPs) rather than phospholipids. For the first time we demonstrate the feasibility of combining vesicular membrane formation and biocatalytic activity with molecular entities of a single class: proteins. The presented self-assembled protein-membrane-based compartments (PMBCs) accommodate either an anabolic reaction, based on free DNA ligase as an example of information transformation processes, or a catabolic process. We present a catabolic process based on a single molecular entity combining an amphiphilic protein with tobacco etch virus (TEV) protease as part of the enclosure of a reaction space and facilitating selective catalytic transformations. Combining compartmentalization and biocatalytic activity by utilizing an amphiphilic molecular building block with and without enzyme functionalization enables new strategies in bottom-up synthetic biology, regenerative medicine, pharmaceutical science and biotechnology.

Comprehensive benchmarking of Markov chain Monte Carlo methods for dynamical systems

BACKGROUND

In quantitative biology, mathematical models are used to describe and analyze biological processes. The parameters of these models are usually unknown and need to be estimated from experimental data using statistical methods. In particular, Markov chain Monte Carlo (MCMC) methods have become increasingly popular as they allow for a rigorous analysis of parameter and prediction uncertainties without the need for assuming parameter identifiability or removing non-identifiable parameters. A broad spectrum of MCMC algorithms have been proposed, including single- and multi-chain approaches. However, selecting and tuning sampling algorithms suited for a given problem remains challenging and a comprehensive comparison of different methods is so far not available.

RESULTS

We present the results of a thorough benchmarking of state-of-the-art single- and multi-chain sampling methods, including Adaptive Metropolis, Delayed Rejection Adaptive Metropolis, Metropolis adjusted Langevin algorithm, Parallel Tempering and Parallel Hierarchical Sampling. Different initialization and adaptation schemes are considered. To ensure a comprehensive and fair comparison, we consider problems with a range of features such as bifurcations, periodical orbits, multistability of steady-state solutions and chaotic regimes. These problem properties give rise to various posterior distributions including uni- and multi-modal distributions and non-normally distributed mode tails. For an objective comparison, we developed a pipeline for the semi-automatic comparison of sampling results.

CONCLUSION

The comparison of MCMC algorithms, initialization and adaptation schemes revealed that overall multi-chain algorithms perform better than single-chain algorithms. In some cases this performance can be further increased by using a preceding multi-start local optimization scheme. These results can inform the selection of sampling methods and the benchmark collection can serve for the evaluation of new algorithms. Furthermore, our results confirm the need to address exploration quality of MCMC chains before applying the commonly used quality measure of effective sample size to prevent false analysis conclusions.

Phase space trajectories and Lyapunov exponents in the dynamics of an alpha-helical protein lattice with intra- and inter-spine interactions

The nonlinear dynamics of intra- and inter-spine interaction models of alpha-helical proteins is investigated by proposing a Hamiltonian using the first quantized operators. Hamilton's equations of motion are derived, and the dynamics is studied by constructing the trajectories and phase space plots in both cases. The phase space plots display a chaotic behaviour in the dynamics, which opens questions about the relationship between the chaos and exciton-exciton and exciton-phonon interactions. This is verified by plotting the Lyapunov characteristic exponent curves.

The study of unfoldable self-avoiding walks — Application to protein structure prediction software

Self-avoiding walks (SAWs) are the source of very difficult problems in probability and enumerative combinatorics. They are of great interest as, for example, they are the basis of protein structure prediction (PSP) in bioinformatics. The authors of this paper have previously shown that, depending on the prediction algorithm, the sets of obtained walk conformations differ: For example, all the SAWs can be generated using stretching-based algorithms whereas only the unfoldable SAWs can be obtained with methods that iteratively fold the straight line. A deeper study of (non-)unfoldable SAWs is presented in this paper. The contribution is first a survey of what is currently known about these sets. In particular, we provide clear definitions of various subsets of SAWs related to pivot moves (unfoldable and non-unfoldable SAWs, etc.) and the first results that we have obtained, theoretically or computationally, on these sets. Then a new theorem on the number of non-unfoldable SAWs is demonstrated. Finally, a list of open questions is provided and the consequences on the PSP problem is proposed.

New ways to boost molecular dynamics simulations

We describe a set of algorithms that allow to simulate dihydrofolate reductase (DHFR, a common benchmark) with the AMBER all-atom force field at 160 nanoseconds/day on a single Intel Core i7 5960X CPU (no graphics processing unit (GPU), 23,786 atoms, particle mesh Ewald (PME), 8.0 Å cutoff, correct atom masses, reproducible trajectory, CPU with 3.6 GHz, no turbo boost, 8 AVX registers). The new features include a mixed multiple time-step algorithm (reaching 5 fs), a tuned version of LINCS to constrain bond angles, the fusion of pair list creation and force calculation, pressure coupling with a "densostat," and exploitation of new CPU instruction sets like AVX2. The impact of Intel's new transactional memory, atomic instructions, and sloppy pair lists is also analyzed. The algorithms map well to GPUs and can automatically handle most Protein Data Bank (PDB) files including ligands. An implementation is available as part of the YASARA molecular modeling and simulation program from www.YASARA.org.

Molecular nonlinear dynamics and protein thermal uncertainty quantification

This work introduces molecular nonlinear dynamics (MND) as a new approach for describing protein folding and aggregation. By using a mode system, we show that the MND of disordered proteins is chaotic while that of folded proteins exhibits intrinsically low dimensional manifolds (ILDMs). The stability of ILDMs is found to strongly correlate with protein energies. We propose a novel method for protein thermal uncertainty quantification based on persistently invariant ILDMs. Extensive comparison with experimental data and the state-of-the-art methods in the field validate the proposed new method for protein B-factor prediction.

IS PROTEIN FOLDING PROBLEM REALLY A NP-COMPLETE ONE? FIRST INVESTIGATIONS

To determine the 3D conformation of proteins is a necessity to understand their functions or interactions with other molecules. It is commonly admitted that, when proteins fold from their primary linear structures to their final 3D conformations, they tend to choose the ones that minimize their free energy. To find the 3D conformation of a protein knowing its amino acid sequence, bioinformaticians use various models of different resolutions and artificial intelligence tools, as the protein folding prediction problem is a NP complete one. More precisely, to determine the backbone structure of the protein using the low resolution models (2D HP square and 3D HP cubic), by finding the conformation that minimizes free energy, is intractable exactly. Both proofs of NP-completeness and the 2D prediction consider that acceptable conformations have to satisfy a self-avoiding walk (SAW) requirement, as two different amino acids cannot occupy a same position in the lattice. It is shown in this document that the SAW requirement considered when proving NP-completeness is different from the SAW requirement used in various prediction programs, and that they are different from the real biological requirement. Indeed, the proof of NP completeness and the predictions in silico consider conformations that are not possible in practice. Consequences of this fact are investigated in this research work.

Molecular modeling of nucleic acid structure: energy and sampling

An overview of computer simulation techniques as applied to nucleic acid systems is presented. This unit discusses methods used to treat energy and to sample representative configurations. Emphasis is placed on molecular mechanics and empirical force fields.

Molecular dynamics simulations of aqueous ions at the liquid-vapor interface accelerated using graphics processors

Molecular dynamics (MD) simulations are a vital tool in chemical research, as they are able to provide an atomistic view of chemical systems and processes that is not obtainable through experiment. However, large-scale MD simulations require access to multicore clusters or supercomputers that are not always available to all researchers. Recently, scientists have returned to exploring the power of graphics processing units (GPUs) for various applications, such as MD, enabled by the recent advances in hardware and integrated programming interfaces such as NVIDIA's CUDA platform. One area of particular interest within the context of chemical applications is that of aqueous interfaces, the salt solutions of which have found application as model systems for studying atmospheric process as well as physical behaviors such as the Hoffmeister effect. Here, we present results of GPU-accelerated simulations of the liquid-vapor interface of aqueous sodium iodide solutions. Analysis of various properties, such as density and surface tension, demonstrates that our model is consistent with previous studies of similar systems. In particular, we find that the current combination of water and ion force fields coupled with the ability to simulate surfaces of differing area enabled by GPU hardware is able to reproduce the experimental trend of increasing salt solution surface tension relative to pure water. In terms of performance, our GPU implementation performs equivalent to CHARMM running on 21 CPUs. Finally, we address possible issues with the accuracy of MD simulaions caused by nonstandard single-precision arithmetic implemented on current GPUs.

Bad Seeds Sprout Perilous Dynamics: Stochastic Thermostat Induced Trajectory Synchronization in Biomolecules

Molecular dynamics simulations starting from different initial conditions are commonly used to mimic the behavior of an experimental ensemble. We show in this article that when a Langevin thermostat is used to maintain constant temperature during such simulations, extreme care must be taken when choosing the random number seeds to prevent statistical correlation among the MD trajectories. While recent studies have shown that stochastically thermostatted trajectories evolving within a single potential basin with identical random number seeds tend to synchronize, we show that there is a synchronization effect even for complex, biologically relevant systems. We demonstrate this effect in simulations of alanine trimer and pentamer and in a simulation of a temperature-jump experiment for peptide folding of a 14-residue peptide. Even in replica-exchange simulations, in which the trajectories are at different temperatures, we find partial synchronization occurring when the same random number seed is employed. We explain this by extending the recent derivation of the synchronization effect for two trajectories in a harmonic well to the case in which the trajectories are at two different temperatures. Our results suggest several ways in which mishandling selection of a pseudorandom number generator initial seed can lead to corruption of simulation data. Simulators can fall into this trap in simple situations such as neglecting to specifically indicate different random seeds in either parallel or sequential restart simulations, utilizing a simulation package with a weak pseudorandom number generator, or using an advanced simulation algorithm that has not been programmed to distribute initial seeds.

Hemoglobin senses body temperature

When aspirating human red blood cells (RBCs) into 1.3 mum pipettes (DeltaP = -2.3 kPa), a transition from blocking the pipette below a critical temperature T(c) = 36.3 +/- 0.3 degrees C to passing it above the T(c) occurred (micropipette passage transition). With a 1.1 mum pipette no passage was seen which enabled RBC volume measurements also above T(c). With increasing temperature RBCs lost volume significantly faster below than above a T(c) = 36.4 +/- 0.7 (volume transition). Colloid osmotic pressure (COP) measurements of RBCs in autologous plasma (25 degrees C < or = T < or = 39.5 degrees C) showed a T (c) at 37.1 +/- 0.2 degrees C above which the COP rapidly decreased (COP transition). In NMR T(1)-relaxation time measurements, the T(1) of RBCs in autologous plasma changed from a linear (r = 0.99) increment below T(c) = 37 +/- 1 degrees C at a rate of 0.023 s/K into zero slope above T(c) (RBC T(1) transition).

IN CONCLUSION

An amorphous hemoglobin-water gel formed in the spherical trail, the residual partial sphere of the aspirated RBC. At T(c), a sudden fluidization of the gel occurs. All changes mentioned above happen at a distinct T(c) close to body temperature. The T(c) is moved +0.8 degrees C to higher temperatures when a D(2)O buffer is used. We suggest a mechanism similar to a "glass transition" or a "colloidal phase transition". At T(c), the stabilizing Hb bound water molecules reach a threshold number enabling a partial Hb unfolding. Thus, Hb senses body temperature which must be inscribed in the primary structure of hemoglobin and possibly other proteins.

Statistical prediction and molecular dynamics simulation

We describe a statistical approach to the validation and improvement of molecular dynamics simulations of macromolecules. We emphasize the use of molecular dynamics simulations to calculate thermodynamic quantities that may be compared to experimental measurements, and the use of a common set of energetic parameters across multiple distinct molecules. We briefly review relevant results from the theory of stochastic processes and discuss the monitoring of convergence to equilibrium, the obtaining of confidence intervals for summary statistics corresponding to measured quantities, and an approach to validation and improvement of simulations based on out-of-sample prediction. We apply these methods to replica exchange molecular dynamics simulations of a set of eight helical peptides under the AMBER potential using implicit solvent. We evaluate the ability of these simulations to quantitatively reproduce experimental helicity measurements obtained by circular dichroism. In addition, we introduce notions of statistical predictive estimation for force-field parameter refinement. We perform a sensitivity analysis to identify key parameters of the potential, and introduce Bayesian updating of these parameters. We demonstrate the effect of parameter updating applied to the internal dielectric constant parameter on the out-of-sample prediction accuracy as measured by cross-validation.

The effect of numerical error on the reproducibility of molecular geometry optimizations

Geometry optimization is one of the most often applied techniques in computational drug discovery. Although geometry optimization routines are generally deterministic, the minimization trajectories can be extremely sensitive to initial conditions, especially in case of larger systems such as proteins. Simple manipulations such as coordinate transformations (translations and rotations), file saving and retrieving, and hydrogen addition can introduce small variations ( approximately 0.001 A) in the starting coordinates which can drastically affect the minimization trajectory. With large systems, optimized geometry differences of up to 1 A RMSD and final energy differences of several kcal/mol can be observed when using many commercially available software packages. Differences in computer platforms can also lead to differences in minimization trajectories. Here we demonstrate how routine structure manipulations can introduce small variations in atomic coordinates, which upon geometry optimization, can give rise to unexpectedly large differences in optimized geometries and final energies. We also show how the same minimizations run on different computer platforms can also lead to different results. The implications of these findings on routine computational chemistry procedures are discussed.

Mechanical response and conformational changes of alpha-actinin domains during unfolding: a molecular dynamics study

Alpha-actinin is a cytoskeleton-binding protein involved in the assembly and regulation of the actin filaments. In this work molecular dynamics method was applied to investigate the mechanical behaviour of the human skeletal muscle alpha-actinin. Five configurations were unfolded at an elongation speed of 0.1 nm/ps in order to investigate the conformational changes occurring during the extension process. Moreover, a sensitivity analysis at different velocities was performed for one of the R2-R3 spectrin-like repeat configuration extracted in order to evaluate the effect of the pulling speed on the mechanical behaviour of the molecule. Two different behaviours were recognized with respect to the pulling speed. In particular, at speed higher than 0.025 nm/ps a continuous rearrangement without evident force peaks was obtained, on the contrary at lower speed evident peaks in the range 500-750 pN were detected. R3 repeat resulted more stable than R2 during mechanical unfolding, due to the lower hydrophobic surface available to the solvent. The characterization of the R2-R3 units can be useful for the development of cytoskeleton network models based on stiffness values obtained by analyses performed at the molecular level.

Hamiltonian dynamics of homopolymer chain models

The Hamiltonian dynamics of chains of nonlinearly coupled particles is numerically investigated in two and three dimensions. Simple, off-lattice homopolymer models are used to represent the interparticle potentials. Time averages of observables numerically computed along dynamical trajectories are found to reproduce results given by the statistical mechanics of homopolymer models. The dynamical treatment, however, indicates a nontrivial transition between regimes of slow and fast phase space mixing. Such a transition is inaccessible to a statistical mechanical treatment and reflects a bimodality in the relaxation of time averages to corresponding ensemble averages. It is also found that a change in the energy dependence of the largest Lyapunov exponent indicates the Theta transition between filamentary and globular polymer configurations, clearly detecting the transition even for a finite number of particles.

Conformational Substates Modulate Hydride Transfer in Dihydrofolate Reductase

In earlier studies of the hydride-transfer reaction catalyzed by dihydrofolate reductase (DHFR) we identified features of the protein correlated with variations in the reaction barrier. We extend the scope of those studies by carrying out potential of mean force (PMF) simulations to determine the hydride-transfer barrier in the wild-type protein as well as the G121V and G121S mutants. While our prior studies focused on the reactant state, our current work addresses the full reaction pathway and directly probes the reactive event. The free energy barriers and structural ensembles resulting from these PMF calculations exhibit the same trends reported in our previous work. Fluctuations present in these simulations also exhibit trends associated with differences in the hydride-transfer barrier height. Moreover, vibrational modes anticipated to promote hydride transfer exhibit larger amplitudes in simulations that generate lowered barriers. The results of our study indicate that discrete basins (substates) on a potential energy landscape of the enzyme give rise to distinct hydride-transfer barriers. We suggest that the long-range effects of mutations at position 121 within DHFR are mediated by differentially preorganized protein environments in the context of distinct substate distributions, with concomitant changes to the dynamic properties of the enzyme.

An emerging concept of biomolecular dynamics and function: applications of NMR & MRI

A new concept of protein dynamics has emerged quite recently, and a crucial link between protein dynamics and function has been largely established using recent NMR techniques in the solution state. Protein structure is governed by the thermodynamic principle and may not necessarily be unique in the solution state. Enzyme catalysis, protein folding or allosteric transition occurs on the microsecond to millisecond time scale, implying that in order to prepare the specific nuclear coordinate for the electronic state transition, a protein must rearrange its nuclear coordinates substantially, and this process may generally take a long period of time almost comparable to that of protein folding. Protein coordinates optimized for the electronic reaction may form an energy state--which may be called an "excited state"--that is thermodynamically distinct from the native state. In contrast, the native state is called a "ground state." Relevant NMR techniques developed recently may also have useful application to MRI, since the critical time scale of various reactions in a living system is also around micro- to milliseconds.

The investigation of the effects of counterions in protein dynamics simulations

Molecular simulations able to exactly represent solvated charged proteins are helpful in understanding protein dynamics, structure and function. In the present study we have used two different starting structures of papain (a typical, stable, globular protein of intermediate net charge) and different modeling procedures to evaluate some effects of counterions in simulations. A number of configurations have been generated and relaxed for each system by various combinations of constrained simulated annealing and molecular dynamics procedures, using the AMBER force field. The analysis of trajectories shows that the simulations of solvated proteins are moderately sensitive to the presence of counterions. However, this sensitivity is highly dependent on the starting model and different procedures of equilibration used. The neutralized systems tend to evince smaller root mean square deviations regardless of the system investigated and the simulation procedure used. The results of parameterized fitting of the simulated structures to the crystallographic data, giving quantitative measure of the total charge influence on the stability of various elements of the secondary structure, revealed a clear scatter of different reactions of various systems' secondary structures to counterions addition: some systems apparently were stabilized when neutralized, while the others were not. Thus, one cannot unequivocally state, despite consideration of specific simulation conditions, whether protein secondary structures are more stable when they have neutralized charges. This suggests that caution should be taken when claiming the stabilizing effect of counterions in simulations other than those involving small, unstable polypeptides or highly charged proteins.

Solvent entropy-driven searching for protein modeling examined and tested in simplified models

Solvent entropy is a force to consider in protein folding and protein design but is difficult to model. It is investigated here in the context of the hp model: Two types of residues, hydrophobic and hydrophilic, are modeled on a lattice. Nine chains and two- and three-dimensional simulations are compared. We show that considering solvent entropy alone, efficient folding of lattice chains (identification of the native fold) can be achieved by an entropy-driven simulation on its own. Moreover, in a detailed comparison over a wide range of parameters, entropy-guided searching outperforms an energy-driven search in the model. The combination of energy- and entropy-driven search yields the most efficient searching. It is compared in detail with the above results, indicating also how this solvent shell model may advantageously be implemented in more complex protein modeling simulations.

Molecular modeling of nucleic acid structure: energy and sampling

An overview of computer simulation techniques as applied to nucleic acid systems is presented. This unit expands an accompanying overview unit (UNIT 7.5) by discussing methods used to treat the energy and sample representative configurations. Emphasis is placed on molecular mechanics and empirical force fields.

Molecular dynamics of the anticodon domain of yeast tRNA(Phe): codon-anticodon interaction

We have studied the effect of codon-anticodon interaction on the structure and dynamics of transfer RNAs using molecular dynamics simulations over a nanosecond time scale. From our molecular dynamical investigations of the solvated anticodon domain of yeast tRNA(Phe) in the presence and absence of the codon trinucleotides UUC and UUU, we find that, although at a gross level the structures are quite similar for the free and the bound domains, there are small but distinct differences in certain parts of the molecule, notably near the Y37 base. Comparison of the dynamics in terms of interatomic or inter-residual distance fluctuation for the free and the bound domains showed regions of enhanced rigidity in the loop region in the presence of codons. Because fluorescence experiments suggested the existence of multiple conformers of the anticodon domain, which interconvert on a much larger time scale than our simulations, we probed the conformational space using five independent trajectories of 500 ps duration. A generalized ergodic measure analysis of the trajectories revealed that at least for this time scale, all the trajectories populated separate parts of the conformational space, indicating a need for even longer simulations or enhanced sampling of the conformational space to give an unequivocal answer to this question.

Polymer principles of protein calorimetric two-state cooperativity

The experimental calorimetric two-state criterion requires the van't Hoff enthalpy DeltaH(vH) around the folding/unfolding transition midpoint to be equal or very close to the calorimetric enthalpy DeltaH(cal) of the entire transition. We use an analytical model with experimental parameters from chymotrypsin inhibitor 2 to elucidate the relationship among several different van't Hoff enthalpies used in calorimetric analyses. Under reasonable assumptions, the implications of these DeltaH(vH)'s being approximately equal to DeltaH(cal) are equivalent: Enthalpic variations among denatured conformations in real proteins are much narrower than some previous lattice-model estimates, suggesting that the energy landscape theory "folding to glass transition temperature ratio" T(f) /T(g) may exceed 6.0 for real calorimetrically two-state proteins. Several popular three-dimensional lattice protein models, with different numbers of residue types in their alphabets, are found to fall short of the high experimental standard for being calorimetrically two-state. Some models postulate a multiple-conformation native state with substantial pre-denaturational energetic fluctuations well below the unfolding transition temperature, or predict a significant post-denaturational continuous conformational expansion of the denatured ensemble at temperatures well above the transition point, or both. These scenarios either disagree with experiments on protein size and dynamics, or are inconsistent with conventional interpretation of calorimetric data. However, when empirical linear baseline subtractions are employed, the resulting DeltaH(vH)/DeltaH(cal)'s for some models can be increased to values closer to unity, and baseline subtractions are found to correspond roughly to an operational definition of native-state conformational diversity. These results necessitate a re-assessment of theoretical models and experimental interpretations.

Unfolding proteins by external forces and temperature: the importance of topology and energetics

Unfolding of proteins by forced stretching with atomic force microscopy or laser tweezer experiments complements more classical techniques using chemical denaturants or temperature. Forced unfolding is of particular interest for proteins that are under mechanical stress in their biological function. For beta-sandwich proteins (a fibronectin type III and an immunoglobulin domain), both of which appear in the muscle protein titin, the results of stretching simulations show important differences from temperature-induced unfolding, but there are common features that point to the existence of folding cores. Intermediates detected by comparing unfolding with a biasing perturbation and a constant pulling force are not evident in temperature-induced unfolding. For an alpha-helical domain (alpha-spectrin), which forms part of the cytoskeleton, there is little commonality in the pathways from unfolding induced by stretching and temperature. Comparison of the forced unfolding of the two beta-sandwich proteins and two alpha-helical proteins (the alpha-spectrin domain and an acyl-coenzyme A-binding protein) highlights important differences within and between protein classes that are related to the folding topologies and the relative stability of the various structural elements.

A molecular dynamics study of the 41-56 beta-hairpin from B1 domain of protein G

The structural and dynamical behavior of the 41-56 beta-hairpin from the protein G B1 domain (GB1) has been studied at different temperatures using molecular dynamics (MD) simulations in an aqueous environment. The purpose of these simulations is to establish the stability of this hairpin in view of its possible role as a nucleation site for protein folding. The conformation of the peptide in the crystallographic structure of the protein GB1 (native conformation) was lost in all simulations. The new equilibrium conformations are stable for several nanoseconds at 300K (>10 ns), 350 K (>6.5 ns), and even at 450 K (up to 2.5 ns). The new structures have very similar hairpin-like conformations with properties in agreement with available experimental nuclear Overhauser effect (NOE) data. The stability of the structure in the hydrophobic core region during the simulations is consistent with the experimental data and provides further evidence for the role played by hydrophobic interactions in hairpin structures. Essential dynamics analysis shows that the dynamics of the peptide at different temperatures spans basically the same essential subspace. The main equilibrium motions in this subspace involve large fluctuations of the residues in the turn and ends regions. Of the six interchain hydrogen bonds, the inner four remain stable during the simulations. The space spanned by the first two eigenvectors, as sampled at 450 K, includes almost all of the 47 different hairpin structures found in the database. Finally, analysis of the hydration of the 300 K average conformations shows that the hydration sites observed in the native conformation are still well hydrated in the equilibrium MD ensemble.

Study of the stability and unfolding mechanism of BBA1 by molecular dynamics simulations at different temperatures

BBA1 is a designed protein that has only 23 residues. It is the smallest protein without disulfide bridges that has a well-defined tertiary structure in solution. We have performed unfolding molecular dynamics simulations on BBA1 and some of its mutants at 300, 330, 360, and 400 K to study their kinetic stability as well as the unfolding mechanism of BBA1. It was shown that the unfolding simulations can provide insights into the forces that stabilize the protein. Packing, hydrophobic interactions, and a salt bridge between Asp12 and Lys16 were found to be important to the protein's stability. The unfolding of BBA1 goes through two major steps: (1) disruption of the hydrophobic core and (2) unfolding of the helix. The beta-hairpin remains stable in the unfolding because of the high stability of the type II' turn connecting the two beta-strands.

Forced unfolding of fibronectin type 3 modules: an analysis by biased molecular dynamics simulations

Titin, an important constituent of vertebrate muscles, is a protein of the order of a micrometer in length in the folded state. Atomic force microscopy and laser tweezer experiments have been used to stretch titin molecules to more than ten times their folded lengths. To explain the observed relation between force and extension, it has been suggested that the immunoglobulin and fibronectin domains unfold one at a time in an all-or-none fashion. We use molecular dynamics simulations to study the forced unfolding of two different fibronectin type 3 domains (the ninth, 9Fn3, and the tenth, 10Fn3, from human fibronectin) and of their heterodimer of known structure. An external biasing potential on the N to C distance is employed and the protein is treated in the polar hydrogen representation with an implicit solvation model. The latter provides an adiabatic solvent response, which is important for the nanosecond unfolding simulation method used here. A series of simulations is performed for each system to obtain meaningful results. The two different fibronectin domains are shown to unfold in the same way along two possible pathways. These involve the partial separation of the "beta-sandwich", an essential structural element, and the unfolding of the individual sheets in a stepwise fashion. The biasing potential results are confirmed by constant force unfolding simulations. For the two connected domains, there is complete unfolding of one domain (9Fn3) before major unfolding of the second domain (10Fn3). Comparison of different models for the potential energy function demonstrates that the dominant cohesive element in both proteins is due to the attractive van der Waals interactions; electrostatic interactions play a structural role but appear to make only a small contribution to the stabilization of the domains, in agreement with other studies of beta-sheet stability. The unfolding forces found in the simulations are of the order of those observed experimentally, even though the speed of the former is more than six orders of magnitude greater than that used in the latter.

Characterization of flexible molecules in solution: the RGDW peptide

Molecular dynamics simulations with adaptive umbrella sampling of the potential energy are used to study conformations of the adhesion peptide RGDW. The peptide is simulated in a box of explicit water. It results in a combination of room temperature (300 K) simulations, in which conformations dominating the average properties of the system are sampled, with high temperature ( approximately 1000 K) simulations in which free energy barriers separating different local minima are crossed efficiently. The simulations with explicit water are compared to simulations of the isolated peptide using different treatments of the electrostatics, and to published experimental data. There is good agreement for data related to the backbone conformation of the peptide. Some discrepancies are evident for data related to side-chain conformations. Together the simulations and experiments provide a description of the RGDW system that is more detailed and reliable than what can be obtained by either simulations or experiments alone.

The early stage of folding of villin headpiece subdomain observed in a 200-nanosecond fully solvated molecular dynamics simulation

A new approach in implementing classical molecular dynamics simulation for parallel computers has enabled a simulation to be carried out on a protein with explicit representation of water an order of magnitude longer than previously reported and will soon enable such simulations to be carried into the microsecond time range. We have used this approach to study the folding of the villin headpiece subdomain, a 36-residue small protein consisting of three helices, from an unfolded structure to a molten globule state, which has a number of features of the native structure. The time development of the solvation free energy, the radius of gyration, and the mainchain rms difference from the native NMR structure showed that the process can be seen as a 60-nsec "burst" phase followed by a slow "conformational readjustment" phase. We found that the burial of the hydrophobic surface dominated the early phase of the folding process and appeared to be the primary driving force of the reduction in the radius of gyration in that phase.

Characterizing global substates of myoglobin

BACKGROUND

The massive amount of information generated from current molecular dynamics simulations makes the data difficult to analyze efficiently. Principal component analysis has been used for almost a century to detect and characterize data relationships and to reduce the dimensionality for problems in many fields. Here, we present an adaptation of principal component analysis using a partial singular value decomposition (SVD) for investigating both the localized and global motions of macromolecules.

RESULTS

Configuration space projections from the SVD analysis of a variety of myoglobin simulations are used to characterize the dynamics of the protein. This technique reveals new dynamical motifs, which quantify proposed hierarchical structures of conformational substates for proteins and provide a means by which configuration space sampling efficiency may be probed. The SVD clearly shows that solvent effects facilitate transitions between global conformational substates for myoglobin molecular dynamics simulations. Lyapunov exponents calculated from the configuration space divergence of 15 trajectories agree with previous predictions for the chaotic behavior of complex protein systems.

CONCLUSIONS

Configuration space projections provide invaluable information about protein motions that would be extremely difficult to obtain otherwise. While the configuration space for myoglobin is quite large, it does have structure. Our analysis of this structure shows that the protein hops between a number of distinct global conformational states, much like the local behavior observed for an individual residue.