Viktor Hornak

Comparison of multiple Amber force fields and development of improved protein backbone parameters.

The ff94 force field that is commonly associated with the Amber simulation package is one of the most widely used parameter sets for biomolecular simulation. After a decade of extensive use and testing, limitations in this force field, such as over‐stabilization of α‐helices, were reported by us and other researchers. This led to a number of attempts to improve these parameters, resulting in a variety of “Amber” force fields and significant difficulty in determining which should be used for a particular application. We show that several of these continue to suffer from inadequate balance between different secondary structure elements. In addition, the approach used in most of these studies neglected to account for the existence in Amber of two sets of backbone φ/ψ dihedral terms. This led to parameter sets that provide unreasonable conformational preferences for glycine. We report here an effort to improve the φ/ψ dihedral terms in the ff99 energy function. Dihedral term parameters are based on fitting the energies of multiple conformations of glycine and alanine tetrapeptides from high level ab initio quantum mechanical calculations. The new parameters for backbone dihedrals replace those in the existing ff99 force field. This parameter set, which we denote ff99SB, achieves a better balance of secondary structure elements as judged by improved distribution of backbone dihedrals for glycine and alanine with respect to PDB survey data. It also accomplishes improved agreement with published experimental data for conformational preferences of short alanine peptides and better accord with experimental NMR relaxation data of test protein systems. Proteins 2006.

Helix movement is coupled to displacement of the second extracellular loop in rhodopsin activation

The second extracellular loop (EL2) of rhodopsin forms a cap over the binding site of its photoreactive 11-cis retinylidene chromophore. A crucial question has been whether EL2 forms a reversible gate that opens upon activation or acts as a rigid barrier. Distance measurements using solid-state 13C NMR spectroscopy between the retinal chromophore and the β4 strand of EL2 show that the loop is displaced from the retinal binding site upon activation, and there is a rearrangement in the hydrogen-bonding networks connecting EL2 with the extracellular ends of transmembrane helices H4, H5 and H6. NMR measurements further reveal that structural changes in EL2 are coupled to the motion of helix H5 and breaking of the ionic lock that regulates activation. These results provide a comprehensive view of how retinal isomerization triggers helix motion and activation in this prototypical G protein–coupled receptor.

Crystallographic structure of the K intermediate of bacteriorhodopsin: conservation of free energy after photoisomerization of the retinal.

The K state, an early intermediate of the bacteriorhodopsin photocycle, contains the excess free energy used for light-driven proton transport. The energy gain must reside in or near the photoisomerized retinal, but in what form has long been an open question. We produced the K intermediate in bacteriorhodopsin crystals in a photostationary state at 100K, with 40% yield, and determined its X-ray diffraction structure to 1.43 A resolution. In independent refinements of data from four crystals, the changes are confined mainly to the photoisomerized retinal. The retinal is 13-cis,15-anti, as known from vibrational spectroscopy. The C13=C14 bond is rotated nearly fully to cis from the initial trans configuration, but the C14-C15 and C15=NZ bonds are partially counter-rotated. This strained geometry keeps the direction of the Schiff base N-H bond vector roughly in the extracellular direction, but the angle of its hydrogen bond with water 402, that connects it to the anionic Asp85 and Asp212, is not optimal. Weakening of this hydrogen bond may account for many of the reported features of the infrared spectrum of K, and for its photoelectric signal, as well as the deprotonation of the Schiff base later in the cycle. Importantly, although 13-cis, the retinal does not assume the expected bent shape of this configuration. Comparison of the calculated energy of the increased angle of C12-C13=C14, that allows this distortion, with the earlier reported calorimetric measurement of the enthalpy gain of the K state indicates that a significant part of the excess energy is conserved in the bond strain at C13.

Location of the Retinal Chromophore in the Activated State of Rhodopsin

Rhodopsin is a highly specialized G protein-coupled receptor (GPCR) that is activated by the rapid photochemical isomerization of its covalently bound 11-cis-retinal chromophore. Using two-dimensional solid-state NMR spectroscopy, we defined the position of the retinal in the active metarhodopsin II intermediate. Distance constraints were obtained between amino acids in the retinal binding site and specific 13C-labeled sites located on the β-ionone ring, polyene chain, and Schiff base end of the retinal. We show that the retinal C20 methyl group rotates toward the second extracellular loop (EL2), which forms a cap on the retinal binding site in the inactive receptor. Despite the trajectory of the methyl group, we observed an increase in the C20-Gly188 (EL2) distance consistent with an increase in separation between the retinal and EL2 upon activation. NMR distance constraints showed that the β-ionone ring moves to a position between Met207 and Phe208 on transmembrane helix H5. Movement of the ring toward H5 was also reflected in increased separation between the Cϵ carbons of Lys296 (H7) and Met44 (H1) and between Gly121 (H3) and the retinal C18 methyl group. Helix-helix interactions involving the H3-H5 and H4-H5 interfaces were also found to change in the formation of metarhodopsin II reflecting increased retinal-protein interactions in the region of Glu122 (H3) and His211 (H5). We discuss the location of the retinal in metarhodopsin II and its interaction with sequence motifs, which are highly conserved across the pharmaceutically important class A GPCR family, with respect to the mechanism of receptor activation.

Using PC clusters to evaluate the transferability of molecular mechanics force fields for proteins

The transferability of molecular mechanics parameters derived for small model systems to larger biopolymers such as proteins can be difficult to assess. Even for small peptides, molecular dynamics simulations are typically too short to sample structures significantly different than initial conformations, making comparison to experimental data questionable. We employed a PC cluster to generate large numbers of native and non‐native conformations for peptides with experimentally measured structural data, one predominantly helical and the other forming a β‐hairpin. These atomic‐detail sets do not suffer from slow convergence, and can be used to rapidly evaluate important force field properties. In this case a suspected bias toward α‐helical conformations in the ff94 and ff99 force fields distributed with the AMBER package was verified. The sets provide critical feedback not only on force field transferability, but may also predict modifications for improvement. Such predictions were used to modify the ff99 parameter set, and the resulting force field was used to test stability and folding of model peptides. Structural behavior during molecular dynamics with the modified force field is found to be very similar to expectations, suggesting that these basis sets of conformations may themselves have significant transferability among force fields. We continue to improve and expand this data set and plan to make it publicly accessible. The calculations involved in this process are trivially parallel and can be performed using inexpensive personal computers with commodity components.

Light activation of rhodopsin: insights from molecular dynamics simulations guided by solid-state NMR distance restraints

Structural restraints provided by solid-state NMR measurements of the metarhodopsin II intermediate are combined with molecular dynamics simulations to help visualize structural changes in the light activation of rhodopsin. Since the timescale for the formation of the metarhodopsin II intermediate (>1 ms) is beyond that readily accessible by molecular dynamics, we use NMR distance restraints derived from 13C dipolar recoupling measurements to guide the simulations. The simulations yield a working model for how photoisomerization of the 11-cis retinylidene chromophore bound within the interior of rhodopsin is coupled to transmembrane helix motion and receptor activation. The mechanism of activation that emerges is that multiple switches on the extracellular (or intradiscal) side of rhodopsin trigger structural changes that converge to disrupt the ionic lock between helices H3 and H6 on the intracellular side of the receptor.

Generation of accurate protein loop conformations through low‐barrier molecular dynamics

Prediction and refinement of protein loop structures are important and challenging tasks for which no general solution has been found. In addition to the accuracy of scoring functions, the main problems reside in (1) insufficient statistical sampling and (2) crossing energy barriers that impede conformational rearrangements of the loop. We approach these two issues by using “low‐barrier molecular dynamics,” a combination of energy smoothing techniques. To address statistical sampling, locally enhanced sampling (LES) is used to produce multiple copies of the loop, thus improving statistics and reducing energy barriers. We introduce a novel extension of LES that can improve local sampling even further through hierarchical subdivision of copies. Even though LES reduces energy barriers, it cannot provide for crossing infinite barriers, which can be problematic when substantial rearrangement of residues is necessary. To permit this kind of loop residue repacking, a “soft‐core” potential energy function is introduced, so that atomic overlaps are temporarily allowed. We tested this new combined methodology to a loop in anti‐influenza antibody Fab 17/9 (7 residues long) and to another loop in the antiprogesterone antibody DB3 (8 residues). In both cases, starting from random conformations, we were able to locate correct loop structures (including sidechain orientations) with heavy‐atom root‐mean‐square deviation (fit to the nonloop region) of ∼1.1 Å in Fab 17/9 and ∼1.8 Å in DB3. We show that the combination of LES and soft‐core potential substantially improves sampling compared to regular molecular dynamics. Moreover, the sampling improvement obtained with this combined approach is significantly better than that provided by either of the two methods alone. Proteins 2003;51:577–590.

Role of group-conserved residues in the helical core of β2-adrenergic receptor

G protein-coupled receptors (GPCRs) belonging to class A contain several highly conserved (>90%) amino acids in their transmembrane helices. Results of mutational studies of these highly conserved residues suggest a common mechanism for locking GPCRs in an inactive conformation and for their subsequent activation upon ligand binding. Recently, a second set of sites in the transmembrane helices has been identified in which amino acids with small side chains, such as Gly, Ala, Ser, Thr, and Cys, are highly conserved (>90%) when considered as a group. These group-conserved residues have not been recognized as having essential structural or functional roles. To determine the role of group-conserved residues in the β2-adrenergic receptor (β2-AR), amino acid replacements guided by molecular modeling were carried out at key positions in transmembrane helices H2–H4. The most significant changes in receptor expression and activity were observed upon replacement of the amino acids Ser-161 and Ser-165 in H4. Substitution at these sites by larger residues lowered the expression and activity of the receptor but did not affect specific binding to the antagonist ligand dihydroalprenolol. A second site mutation, V114A, rescued the low expression of the S165V mutant. Substitution of other group-conserved residues in H2–H4 by larger amino acids lowered receptor activity in the order Ala-128, Ala-76, Ser-120, and Ala-78. Together these data provide comprehensive analysis of group-conserved residues in a class A GPCR and allow insights into the roles of these residues in GPCR structure and function.

Reconciling the solution and X-ray structures of the villin headpiece helical subdomain: molecular dynamics simulations and double mutant cycles reveal a stabilizing cation-pi interaction.

The 36-residue helical subdomain of the villin headpiece, HP36, is one of the smallest cooperatively folded proteins, folding on the microsecond time scale. The domain is an extraordinarily popular model system for both experimental and computational studies of protein folding. The structure of HP36 has been determined using X-ray crystallography and NMR spectroscopy, with the resulting structures exhibiting differences in helix packing, van der Waals contacts, and hydrogen bonding. It is important to determine the solution structure of HP36 with as much accuracy as possible since this structure is widely used as a reference for simulations and experiments. We complement the existing data by using all-atom molecular dynamics simulations with explicit solvent to evaluate which of the experimental models is the better representation of HP36 in solution. After simulation for 50 ns initiated with the NMR structure, we observed that the protein spontaneously adopts structures with a backbone conformation, core packing, and C-capping motif on the third helix that are more consistent with the crystal structure. We also examined hydrogen bonding and side chain packing interactions between D44 and R55 and between F47 and R55, respectively, which were observed in the crystal structure but not in the NMR-based solution structure. Simulations showed large fluctuations in the distance between D44 and R55, while the distance between F47 and R55 remained stable, suggesting the formation of a cation-pi interaction between those residues. Experimental double mutant cycles confirmed that the F47-R55 pair has a larger energetic coupling than the D44-R55 interaction. Overall, these combined experimental and computational studies show that the X-ray crystal structure is the better reference structure for HP36 in solution at neutral pH. Our analysis also shows how detailed molecular dynamics simulations combined with experimental validation can help bridge the gap between NMR and crystallographic methods.

Structure–Kinetic Relationships of Passive Membrane Permeation from Multiscale Modeling

Passive membrane permeation of small molecules is essential to achieve the required absorption, distribution, metabolism, and excretion (ADME) profiles of drug candidates, in particular intestinal absorption and transport across the blood-brain barrier. Computational investigations of this process typically involve either building QSAR models or performing free energy calculations of the permeation event. Although insightful, these methods rarely bridge the gap between computation and experiment in a quantitative manner, and identifying structural insights to apply toward the design of compounds with improved permeability can be difficult. In this work, we combine molecular dynamics simulations capturing the kinetic steps of permeation at the atomistic level with a dynamic mechanistic model describing permeation at the in vitro level, finding a high level of agreement with experimental permeation measurements. Calculation of the kinetic rate constants determining each step in the permeation event allows derivation of structure-kinetic relationships of permeation. We use these relationships to probe the structural determinants of membrane permeation, finding that the desolvation/loss of hydrogen bonding required to leave the membrane partitioned position controls the membrane flip-flop rate, whereas membrane partitioning determines the rate of leaving the membrane.