Brownian dynamics simulations of protein folding: Access to milliseconds time scale and beyond

Protein folding occurs on a time scale ranging from milliseconds to minutes for a majority of proteins. Computer simulation of protein folding, from a random configuration to the native structure, is nontrivial owing to the large disparity between the simulation and folding time scales. As an effort to overcome this limitation, simple models with idealized protein subdomains, e.g., the diffusion–collision model of Karplus and Weaver, have gained some popularity. We present here new results for the folding of a four-helix bundle within the framework of the diffusion–collision model. Even with such simplifying assumptions, a direct application of standard Brownian dynamics methods would consume 10,000 processor-years on current supercomputers. We circumvent this difficulty by invoking a special Brownian dynamics simulation. The method features the calculation of the mean passage time of an event from the flux overpopulation method and the sampling of events that lead to productive collisions even if their probability is extremely small (because of large free-energy barriers that separate them from the higher probability events). Using these developments, we demonstrate that a coarse-grained model of the four-helix bundle can be simulated in several days on current supercomputers. Furthermore, such simulations yield folding times that are in the range of time scales observed in experiments.

Folding protein models with a simple hydrophobic energy function: The fundamental importance of monomer inside/outside segregation

The present study explores a “hydrophobic” energy function for folding simulations of the protein lattice model. The contribution of each monomer to conformational energy is the product of its “hydrophobicity” and the number of contacts it makes, i.e., E(h⃗, c⃗) = −Σi=1N cihi = −(h⃗.c⃗) is the negative scalar product between two vectors in N-dimensional cartesian space: h⃗ = (h1, … , hN), which represents monomer hydrophobicities and is sequence-dependent; and c⃗ = (c1, … , cN), which represents the number of contacts made by each monomer and is conformation-dependent. A simple theoretical analysis shows that restrictions are imposed concomitantly on both sequences and native structures if the stability criterion for protein-like behavior is to be satisfied. Given a conformation with vector c⃗, the best sequence is a vector h⃗ on the direction upon which the projection of c⃗ − c̄⃗ is maximal, where c̄⃗ is the diagonal vector with components equal to c̄, the average number of contacts per monomer in the unfolded state. Best native conformations are suggested to be not maximally compact, as assumed in many studies, but the ones with largest variance of contacts among its monomers, i.e., with monomers tending to occupy completely buried or completely exposed positions. This inside/outside segregation is reflected on an apolar/polar distribution on the corresponding sequence. Monte Carlo simulations in two dimensions corroborate this general scheme. Sequences targeted to conformations with large contact variances folded cooperatively with thermodynamics of a two-state transition. Sequences targeted to maximally compact conformations, which have lower contact variance, were either found to have degenerate ground state or to fold with much lower cooperativity.

Vibrational population relaxation of carbon monoxide in the heme pocket of photolyzed carbonmonoxy myoglobin: Comparison of time-resolved mid-IR absorbance experiments and molecular dynamics simulations

The vibrational energy relaxation of carbon monoxide in the heme pocket of sperm whale myoglobin was studied by using molecular dynamics simulation and normal mode analysis methods. Molecular dynamics trajectories of solvated myoglobin were run at 300 K for both the δ- and ɛ-tautomers of the distal His-64. Vibrational population relaxation times of 335 ± 115 ps for the δ-tautomer and 640 ± 185 ps for the ɛ-tautomer were estimated by using the Landau–Teller model. Normal mode analysis was used to identify those protein residues that act as the primary “doorway” modes in the vibrational relaxation of the oscillator. Although the CO relaxation rates in both the ɛ- and δ-tautomers are similar in magnitude, the simulations predict that the vibrational relaxation of the CO is faster in the δ-tautomer with the distal His playing an important role in the energy relaxation mechanism. Time-resolved mid-IR absorbance measurements were performed on photolyzed carbonmonoxy hemoglobin (Hb13CO). From these measurements, a T1 time of 600 ± 150 ps was determined. The simulation and experimental estimates are compared and discussed.

Molecular dynamics and free-energy calculations applied to affinity maturation in antibody 48G7

We investigated the relative free energies of hapten binding to the germ line and mature forms of the 48G7 antibody Fab fragments by applying a continuum model to structures sampled from molecular dynamics simulations in explicit solvent. Reasonable absolute and very good relative free energies were obtained. As a result of nine somatic mutations that do not contact the hapten, the affinity-matured antibody binds the hapten >104 tighter than the germ line antibody. Energetic analysis reveals that van der Waals interactions and nonpolar contributions to solvation are similar and drive the formations of both the germ line and mature antibody–hapten complexes. Affinity maturation of the 48G7 antibody therefore appears to occur through reorganization of the combining site geometry in a manner that optimizes the balance of gaining favorable electrostatic interactions with the hapten and losing those with solvent during the binding process. As reflected by lower rms fluctuations in the antibody–hapten complex, the mature complex undergoes more restricted fluctuations than the germ line complex. The dramatically increased affinity of the 48G7 antibody over its germ line precursor is thus made possible by electrostatic optimization.

Folding simulations of a three-stranded antiparallel β-sheet peptide

Protein folding is a grand challenge of the postgenomic era. In this paper, 58 folding events sampled during 47 molecular dynamics trajectories for a total simulation time of more than 4 μs provide an atomic detail picture of the folding of a 20-residue synthetic peptide with a stable three-stranded antiparallel β-sheet fold. The simulations successfully reproduce the NMR solution conformation, irrespective of the starting structure. The sampling of the conformational space is sufficient to determine the free energy surface and localize the minima and transition states. The statistically predominant folding pathway involves the formation of contacts between strands 2 and 3, starting with the side chains close to the turn, followed by association of the N-terminal strand onto the preformed 2–3 β-hairpin. The folding mechanism presented here, formation of a β-hairpin followed by consolidation, is in agreement with a computational study of the free energy surface of another synthetic three-stranded antiparallel β-sheet by Bursulaya and Brooks [(1999) J. Am. Chem. Soc. 121, 9947–9951]. Hence, it might hold in general for antiparallel β-sheets with short turns.

Quantification of the hydrophobic interaction by simulations of the aggregation of small hydrophobic solutes in water

The hydrophobic interaction, the tendency for nonpolar molecules to aggregate in solution, is a major driving force in biology. In a direct approach to the physical basis of the hydrophobic effect, nanosecond molecular dynamics simulations were performed on increasing numbers of hydrocarbon solute molecules in water-filled boxes of different sizes. The intermittent formation of solute clusters gives a free energy that is proportional to the loss in exposed molecular surface area with a constant of proportionality of 45 ± 6 cal/mol⋅Å2. The molecular surface area is the envelope of the solute cluster that is impenetrable by solvent and is somewhat smaller than the more traditional solvent-accessible surface area, which is the area transcribed by the radius of a solvent molecule rolled over the surface of the cluster. When we apply a factor relating molecular surface area to solvent-accessible surface area, we obtain 24 cal/mol⋅Å2. Ours is the first direct calculation, to our knowledge, of the hydrophobic interaction from molecular dynamics simulations; the excellent qualitative and quantitative agreement with experiment proves that simple van der Waals interactions and atomic point-charge electrostatics account for the most important driving force in biology.

Computer simulations of enzyme catalysis: Finding out what has been optimized by evolution

The origin of the catalytic power of enzymes is discussed, paying attention to evolutionary constraints. It is pointed out that enzyme catalysis reflects energy contributions that cannot be determined uniquely by current experimental approaches without augmenting the analysis by computer simulation studies. The use of energy considerations and computer simulations allows one to exclude many of the popular proposals for the way enzymes work. It appears that the standard approaches used by organic chemists to catalyze reactions in solutions are not used by enzymes. This point is illustrated by considering the desolvation hypothesis and showing that it cannot account for a large increase in kcat relative to the corresponding kcage for the reference reaction in a solvent cage. The problems associated with other frequently invoked mechanisms also are outlined. Furthermore, it is pointed out that mutation studies are inconsistent with ground state destabilization mechanisms. After considering factors that were not optimized by evolution, we review computer simulation studies that reproduced the overall catalytic effect of different enzymes. These studies pointed toward electrostatic effects as the most important catalytic contributions. The nature of this electrostatic stabilization mechanism is far from being obvious because the electrostatic interaction between the reacting system and the surrounding area is similar in enzymes and in solution. However, the difference is that enzymes have a preorganized dipolar environment that does not have to pay the reorganization energy for stabilizing the relevant transition states. Apparently, the catalytic power of enzymes is stored in their folding energy in the form of the preorganized polar environment.

Unraveling principles of lead discovery: from unfrustrated energy landscapes to novel molecular anchors.

The search for novel leads is a critical step in the drug discovery process. Computational approaches to identify new lead molecules have focused on discovering complete ligands by evaluating the binding affinity of a large number of candidates, a task of considerable complexity. A new computational method is introduced in this work based on the premise that the primary molecular recognition event in the protein binding site may be accomplished by small core fragments that serve as molecular anchors, providing a structurally stable platform that can be subsequently tailored into complete ligands. To fulfill its role, we show that an effective molecular anchor must meet both the thermodynamic requirement of relative energetic stability of a single binding mode and its consistent kinetic accessibility, which may be measured by the structural consensus of multiple docking simulations. From a large number of candidates, this technique is able to identify known core fragments responsible for primary recognition by the FK506 binding protein (FKBP-12), along with a diverse repertoire of novel molecular cores. By contrast, absolute energetic criteria for selecting molecular anchors are found to be promiscuous. A relationship between a minimum frustration principle of binding energy landscapes and receptor-specific molecular anchors in their role as "recognition nuclei" is established, thereby unraveling a mechanism of lead discovery and providing a practical route to receptor-biased computational combinatorial chemistry.