Search

HIV-1 reverse transcriptase is a critical drug target for HIV treatment, and understanding the exact mechanisms of its function and inhibition would significantly accelerate the development of new anti-HIV drugs. It is well known that structure plays a critical role in protein function, but for reverse transcriptase, structural information has proven to be insufficient – despite enormous effort – to explain the mechanism of inhibition and drug resistance of non-nucleoside reverse transcriptase inhibitors. We hypothesize that the missing link is dynamics, information about the motions of the system. However, many of the techniques that give the best information about dynamics, such as solution NMR and molecular dynamics simulations, cannot be easily applied to a protein as large as reverse transcriptase. As an alternative, we combine elastic network modeling with simultaneous hierarchical clustering of structural and dynamic data. We present an extensive survey of the dynamics of reverse transcriptase bound to a variety of ligands and with a number of mutations, revealing a novel mechanism for drug resistance to non-nucleoside reverse transcriptase inhibitors. Hydrophobic core mutations restore active-state motion to multiple functionally significant regions of HIV-1 RT. This model arises out of a combination of structural and dynamic information, rather than exclusively from one or the other.

The backbone dihedral parameters of the Amber RNA force field were improved by fitting using multiple linear regression to potential energies determined by quantum chemistry calculations. Five backbone and four glycosidic dihedral parameters were fit simultaneously to reproduce the potential energies determined by a high-level density functional theory calculation (B97D3 functional with the AUG-CC-PVTZ basis set). Umbrella sampling was used to determine conformational free energies along the dihedral angles, and these better agree with the population of conformations observed in the protein data bank for the new parameters than for the conventional parameters. Molecular dynamics simulations performed on a set of hairpin loops, duplexes and tetramers with the new parameter set show improved modeling for the structures of tetramers CCCC, CAAU, and GACC, and an RNA internal loop of noncanonical pairs, as compared to the conventional parameters. For the tetramers, the new parameters largely avoid the incorrect intercalated structures that dominate the conformational samples from the conventional parameters. For the internal loop, the major conformation solved by NMR is stable with the new parameters, but not with the conventional parameters. The new force field performs similarly to the conventional parameters for the UUCG and GCAA hairpin loops and the [U(UA)6A]2 duplex.

Understanding the functions of biomolecules requires insight not only from structures but from dynamics as well. Often, the most interesting processes occur on time scales too slow for exploration by conventional molecular dynamics (MD) simulations. For this reason, alternative computational methods such as elastic network models (ENMs) have become increasingly popular. These simple, coarse-grained models represent molecules as beads connected by harmonic springs; the system's motions are solved analytically by normal-mode analysis. In the past few years, many different formalisms for performing ENM calculations have emerged, and several have been optimized using all-atom MD simulations. In contrast to other studies, we have compared the various formalisms in a systematic, quantitative way. In this study, we optimize many ENM functional forms using a uniform data set containing only long (>1 μs) all-atom MD simulations. Our results show that all models once optimized produce spring constants for immediate neighboring residues that are orders of magnitude stiffer than more distal contacts. In addition, the statistical significance of ENM performance varied with model resolution. We also show that fitting long trajectories does not improve ENM performance due to a problem inherent in all network models tested: they underestimate the relative importance of the most concerted motions. Finally, we characterize ENMs' resilience by tessellating the parameter space to show that broad ranges of parameters produce similar quality predictions. Taken together, our data reveal that the choice of spring function and parameters are not vital to the performance of a network model and that simple parameters can by derived "by hand" when no data are available for fitting, thus illustrating the robustness of these models.

We present a detailed analysis of the behavior of the highly flexible post- translational lipid modifications of rhodopsin from multiple-microsecond all-atom molecular dynamics simulations. Rhodopsin was studied in a realistic membrane environment that includes cholesterol, as well as saturated and polyunsaturated lipids with phosphocholine and phosphoethanolamine headgroups. The simulation reveals striking differences between the palmitoylations at Cys322 and Cys323 as well as between the palmitoyl chains and the neighboring lipids. Notably the palmitoyl group at Cys322 shows considerably greater contact with helix H1 of rhodopsin, yielding frequent chain upturns with longer reorientational correlation times, and relatively low order parameters. While the palmitoylation at Cys323 makes fewer protein contacts and has increased order compared to Cys322, it nevertheless exhibits greater flexibility with smaller order parameters than the stearoyl chains of the surrounding lipids. The dynamical structure of the palmitoylationsas well as their extensive fluctuationssuggests a complex function for the post-translational modifications in rhodopsin and potentially other G protein-coupled receptors, going beyond their role as membrane anchoring elements. Rather, we propose that the palmitoylation at Cys323 has a potential role as a lipid anchor, whereas the palmitoyl−protein interaction observed for Cys322 suggests a more specific interaction that affects the stability of the dark state of rhodopsin.

Molecular dynamics (MD) is a powerful tool for understanding the fluctuations of biomolecular systems. It is, however, subject to statistical errors in its sampling of the underlying distribution of states. One must understand these errors in order to draw meaningful conclusions from the simulation. This is becoming ever more critical as MD simulations of even larger systems are attempted. We present here a new method for determining the extent of convergence that relies on measures of the fluctuation space sampled by the simulation without any a priori knowledge of states or partitioning of the configuration space. This method reveals long correlation times, even for simple systems, and suggests caution when interpreting macromolecular simulations. We also compare this method with previous efforts to characterize the sampling of MD simulation.

High-resolution structural information is routinely available for soluble proteins, largely from x-ray crystallography and solution nuclear magnetic resonance (NMR). However, these techniques are far harder to apply to membrane-bound proteins; membrane proteins are generally reluctant to crystallize, and the need for a surrounding lipid matrix generally means that NMR must be performed under solid-state as opposed to solution conditions. The latter imposes some restrictions on the kinds of information that can be readily extracted from experiments, and one often must rely on experimental methods such as residual dipolar coupling and chemical shift anisotropy, which yield information about the orientations of specific moieties relative to the magnetic field. Although extremely valuable, this sort of data can be harder to interpret than the nuclear-Overhauser-effect-based distance restraints one customarily sees with soluble proteins. In this issue of Biophysical Journal, De Simone et al. (1) describe their efforts to build structural models from this kind of data by combining experimentally generated restraints with all-atom molecular-dynamics (MD) simulations in explicit lipid bilayers. Specifically, they implemented restraints based on two kinds of solid-state NMR data—chemical shift anisotropy and amide dipolar coupling—and used these restraints to drive calculations on two membrane protein systems, sarcolipin and phospholamban. In both cases, the calculations produced well-formed, reasonable-looking structures consistent with the experimental restraints (and in one case, other experimental data not explicitly used in the simulations). Although other researchers (particularly Im et al. (2), De Simone et al. (3), and Richter et al. (4)) have worked to incorporate NMR observables as restraints into MD simulations, and others have recently described the theoretical basis justifying such an approach (5), this work has a couple of significant strengths. Most notably, using a formalism similar to one presented previously (2), they contrast the behavior of single-trajectory restraints and ensemble restraints; in the latter case, multiple trajectories are run simultaneously, and only the instantaneous average over all trajectories is restrained. The primary advantage of this approach comes when the experimental signal is the result of diverse structures, as opposed to fluctuations about a single state. In any case, it is important to note that this approach involves comparing—or driving—the simulation by comparison to the explicit NMR observable (e.g., the observed dipolar couplings), rather than to the interpretation of those observables (a particular angle or orientation). Previous work on x-ray scattering (6,7) and solid-state 2H NMR (8,9) has shown that taking the former approach is crucial, because where experimentalists are often forced to make sometimes oversimplified assumptions about the conformational ensemble underlying the signals to interpret them, MD simulations explicitly sample those fluctuations. There are two major challenges when trying to incorporate experimental NMR data into a simulation. The first is one of degeneracy. The experimental observables used in this kind of approach typically depend on the sin2 or cos2 of the angle made by some vector (e.g., a bond vector) with the magnetic field, meaning that for any given observable, there are four possible orientations that would produce the same value. This in principle requires a group-theoretic approach to enumerate all possible solutions (10), and if the initial structure used to seed the calculation is sufficiently far from the correct one, the wrong minimum may be sampled. The second challenge is the mismatch between the timescale and time-resolution of solid-state NMR and MD. Every solid-state NMR experiment simultaneously measures two kinds of averages: ensemble and temporal. The former is relatively straightforward: NMR measures the response of the totality of the sample contents, and as such reflects the average over all of the relevant moieties in the system. Even in the case of very narrow selective labeling experiments (10), this means that many millions of particles are simultaneously averaged. By contrast, a MD simulation may have only one (in the case of a membrane protein experiment) or a few hundred (if looking at lipids) signals to average over. However, because these effects are generally additive, one can usually rely on the ergodic hypothesis (the statement that in the limit of infinite sampling, a time average and an ensemble average are equivalent) to rescue the comparison. Although achieving ergodic sampling is a major challenge (11,12), improvements in computer hardware and sampling techniques are making this requirement more and more reasonable. However, NMR also implicitly performs time averaging, due to the finite shutter-speed of the method. We are all familiar with the blurring that occurs when a photographic subject moves while the shutter is open; artistic photographers can take advantage of the effect to produce striking images, but the analogous process is problematic for experimentalists. In the case of NMR, there are several physical origins of this averaging, including the finite duration of the various radio-frequency pulses and mixing times, as well as the nature of the specific phenomena measured (e.g., magnetization transfer, where the rate is related to a correlation function). By contrast, the time-resolution of a MD simulation is specified by the time step (generally 2 fs for all-atom calculations), with no averaging at all. At any point in time (specified with infinite precision), all information about the state of the system (all particle positions and velocities) is fully known. Moreover, experimentally based restraints are generally derived directly from their conformational dependence, without accounting for this averaging. This means that if NMR data is used to drive a MD simulation, there is a risk of a temporal mismatch between the two that could produce odd results. There are three distinct cases that determine what sort of effects this kind of averaging can cause, shown pictorially in Fig. 1. Figure 1 Timescales for exchange and the resulting signals. Three cases are demonstrated: a single state, a pair of states undergoing slow exchange, and a pair of states undergoing fast exchange. In the first case, there is in essence only one state (at least as far as the observable in question is concerned); in this case, the minimum energy of the restraint will correspond to the most populated state, and the system will fluctuate about it realistically. In this case, the simple single-system protocol discussed by De Simone et al. (1) will likely be successful. The second case posits the existence of two distinct states, with different signals, that exchange slowly relative to the NMR timescale. In this case, the observable will also have two distinct peaks, such that in all likelihood there is no single structure consistent with all of the data. For example, an ensemble of unrestrained simulations of rhodopsin showed that the solid-state NMR spectra used to measure the orientation of its ligand were best explained by an ensemble of structures, as opposed to a single state (8,9,13). In this case, the single restrained trajectory approach will fail to reproduce the data. However, an ensemble-based approach, such as that proposed here by De Simone et al. (1) and elsewhere by Im et al. (2), is needed. Here, instead of running a single trajectory, a set of independent trajectories is run simultaneously, and the restraints are applied to the ensemble as a whole as opposed to each trajectory independently. As a result, the ensemble is far more likely to be able to capture the effects of multiple states on the observable spectra. Even here, care must be taken to ensure diverse starting structures, or all of the trajectories may initially fall into the same free energy well; while the restraints should eventually cause some of the trajectories to find the second state, such evolution could be slow if the barriers between the states are high. However, in the third case, where the two states exchange rapidly, the situation is more complex. In this case, the signals for the two states will not appear as two distinct peaks, but rather as a single peak in between. In this case (illustrated in the bottom panel of Fig. 1), naively applied restraints would push the system toward the center, which should not be populated significantly. The ensemble-based problem will not fix this problem, because the issue is one of averaging, not sampling; the restraint is applied instantaneously, but the signal on which it is based is itself an average. To handle this case, a new strategy will likely be necessary. On the one hand, one could refine the experiments to resolve the two states despite their rapid exchange, at which point the ensemble-based method would again be promising. Alternatively, one could explicitly incorporate time-averaging into the restraint, perhaps by applying the restraints to a running average of structures from the molecular dynamics; one could imagine that the formalism would resemble self-guided molecular dynamics (14). Despite this reservation, it is clear that the work of De Simone et al. (1) is of significant value. MD simulations and solid-state NMR are two of the most powerful techniques for exploring the structure and dynamics of membrane proteins, so methods to couple them have the potential for great synergy.

Low-conductance, highly calcium-selective channels formed by the Orai proteins exist as store-operated CRAC channels and store-independent, arachidonic acid–activated ARC channels. Both are activated by stromal interaction molecule 1 (STIM1), but CRAC channels are activated by STIM1 located in the endoplasmic reticulum membrane, whereas ARC channels are activated by the minor plasma membrane–associated pool of STIM1. Critically, maximally activated CRAC channel and ARC channel currents are completely additive within the same cell, and their selective activation results in their ability to each induce distinct cellular responses. We have previously shown that specific ARC channel activation requires a PKA-mediated phosphorylation of a single threonine residue (Thr(389)) within the cytoplasmic region of STIM1. Here, examination of the molecular basis of this phosphorylation-dependent activation revealed that phosphorylation of the Thr(389) residue induces a significant structural change in the STIM1-Orai–activating region (SOAR) that interacts with the Orai proteins, and it is this change that determines the selective activation of the store-independent ARC channels versus the store-operated CRAC channels. In conclusion, our findings reveal the structural changes underlying the selective activation of STIM1-induced CRAC or ARC channels that determine the specific stimulation of these two functionally distinct Ca(2+) entry pathways.

Fengycin is a class of antimicrobial fungicide, synthesized by the bacterial family, Bacillus, which functions by damaging the fungus’ cell membrane. Thus it is cumbersome for the bacteria to develop resistance against it. This makes fengycin a potential drug candidate and comprehending its mechanism of action is a crucial step in drug development. Previously, it has been observed in coarse-grained simulations that fengycins form clusters in model fungal membranes (POPC) but not in model bacterial membranes (POPE:POPG). Our hypothesis is that formation of clusters plays a crucial step for fengycin disrupting the fungal membranes. However, most recent all-atom simulations results suggest that preference for the aggregation process in one membrane over the other is not statistically significant, although this may be more due to slow convergence than a true lack of selectivity. Accordingly, we are using enhanced sampling techniques are being used to compute the free energy change during aggregation in the two membrane types. In addition, it is important to study the effects of sterols on fengycin action, since one of the major differences between mammalian and fungal membranes is the presence of cholesterol versus ergosterol in the outer membrane. Therefore, coarse-grained simulations of lipids with cholesterol will show how these sterols can inhibit or promote the action of fengycin. These results will lead us to a more coherent understanding of the function of this unique biological molecule.

We have developed LOOS (Lightweight Object Oriented Structure-analysis) as a tool for making new tools to analyze molecular simulations. LOOS is an object-oriented library designed to facilitate the rapid development of new methods for structural analysis. LOOS is written in C++ and is easily extensible, only requiring knowledge of 4 core classes. A Python interface is also available, further facilitating rapid development of analysis tools and broadening the LOOS community. LOOS reads the native file formats of most common simulation packages and can write NAMD formats (PDB and DCD) and Gromacs XTC. A dynamic atom selection language, based on C expression syntax, is included and is easily accessible via a single function call. LOOS includes over 140 pre-built tools for common structural analysis tasks, including analyzing simulation convergence, 3D histograms, and elastic network models. Python-based packages include the optimal membrane generator and tools for the Voronoi tesselation of membranes. In addition, tool templates for common analysis work-flows are included to facilitate the development of new tools. LOOS is available for download at http://loos.sourceforge.net and is distributed under the GPLv3 license.

Rhodopsin is currently the only available atomic-resolution template for understanding biological functions of the G protein-coupled receptor (GPCR) family. The structural basis for the phenomenal dark state stability of 11-cis-retinal bound to rhodopsin and its ultrafast photoreaction are active topics of research. In particular, the beta-ionone ring of the retinylidene inverse agonist is crucial for the activation mechanism. We analyzed a total of 23 independent, 100 ns all-atom molecular dynamics simulations of rhodopsin embedded in a lipid bilayer in the microcanonical (N,V,E) ensemble. Analysis of intramolecular fluctuations predicts hydrogen-out-of-plane (HOOP) wagging modes of retinal consistent with those found in Raman vibrational spectroscopy. We show that sampling and ergodicity of the ensemble of simulations are crucial for determining the distribution of conformers of retinal bound to rhodopsin. The polyene chain is rigidly locked into a single, twisted conformation, consistent with the function of retinal as an inverse agonist in the dark state. Most surprisingly, the beta-ionone ring is mobile within its binding pocket; interactions are non-specific and the cavity is sufficiently large to enable structural heterogeneity. We find that retinal occupies two distinct conformations in the dark state, contrary to most previous assumptions. The beta-ionone ring can rotate relative to the polyene chain, thereby populating both positively and negatively twisted 6-s-cis enantiomers. This result, while unexpected, strongly agrees with experimental solid-state (2)H NMR spectra. Correlation analysis identifies the residues most critical to controlling mobility of retinal; we find that Trp265 moves away from the ionone ring prior to any conformational transition. Our findings reinforce how molecular dynamics simulations can challenge conventional assumptions for interpreting experimental data, especially where existing models neglect conformational fluctuations.

The central question in evaluating almost any result from a molecular dynamics simulation is whether the calculation has converged. Unfortunately, assessing the ergodicity of a single trajectory is very difficult to do. In this work, we assess the sampling of molecular dynamics simulations of the membrane protein rhodopsin by comparing the results from 26 independent trajectories, each run for 100 ns. By examining principal components and cluster populations, we show that rhodopsin's fluctuations are not well described by 100 ns of dynamics, and that the sampling is not fully converged even for individual loops. The results serve as a reminder of the caution required when interpreting molecular dynamics simulations of macromolecules.

Evidence from a variety of spectroscopic probes indicates that (phi, psi) values corresponding to the left-handed polyproline II helix (P(II)) are preferred for short alanine-based peptides in water. On the basis of results from theoretical studies, it is believed that the observed preference is dictated by favorable peptide-solvent interactions, which are realized through formation of optimal hydrogen-bonding water bridges between peptide donor and acceptor groups. In the present study, we address this issue explicitly by analyzing the hydration structure and thermodynamics of 16 low-energy conformers of the alanine dipeptide (N-acetylalanine-N'-methylamide) in liquid water. Monte Carlo simulations in the canonical ensemble were performed under ambient conditions with all-atom OPLS parameters for the alanine dipeptide and the TIP5P model for water. We find that the number of hydrogen-bonded water molecules connecting the peptide group donor and acceptor atoms has no effect on the solvation thermodynamics. Instead, the latter are determined by the work done to fully hydrate the peptide. This work is minimal for conformations that are characterized by a minimal overlap of the primary hydration shells around the peptide donor and acceptor atoms. As a result, peptide-solvent interactions favor "compact" conformations that do not include P(II)-like geometries. Our main conclusion is that the experimentally observed preference for P(II) does not arise due to favorable direct interactions between the peptide and water molecules. Instead, the latter act to unmask underlying conformational preferences that are a consequence of minimizing intrapeptide steric conflicts.

Thermodynamic measurements of the solvation of salts and electrolytes are relatively straightforward, but it is not possible to separate total solvation free energies into distinct cation and anion contributions without reference to an additional extrathermodynamic assumption. The present work attempts to resolve this difficulty using molecular dynamics simulations with the AMOEBA polarizable force field and perturbation techniques to directly compute absolute solvation free energies for potassium, sodium, and chloride ions in liquid water and formamide. Corresponding calculations are also performed with two widely used nonpolarizable force fields. The simulations with the polarizable force field accurately reproduce in vacuo quantum mechanical results, experimental ion-cluster solvation enthalpies, and experimental solvation free energies for whole salts, while the other force fields do not. The results indicate that calculations with a polarizable force field can capture the thermodynamics of ion solvation and that the solvation free energies of the individual ions differ by several kilocalories from commonly cited values.

In this study, we applied a comprehensive G protein-coupled receptor-Gαi protein chemical cross-linking strategy to map the cannabinoid receptor subtype 2 (CB2)-Gαi interface and then used molecular dynamics simulations to explore the dynamics of complex formation. Three cross-link sites were identified using LC-MS/MS and electrospray ionization-MS/MS as follows: 1) a sulfhydryl cross-link between C3.53(134) in TMH3 and the Gαi C-terminal i-3 residue Cys-351; 2) a lysine cross-link between K6.35(245) in TMH6 and the Gαi C-terminal i-5 residue, Lys-349; and 3) a lysine cross-link between K5.64(215) in TMH5 and the Gαi α4β6 loop residue, Lys-317. To investigate the dynamics and nature of the conformational changes involved in CB2·Gi complex formation, we carried out microsecond-time scale molecular dynamics simulations of the CB2 R*·Gαi1β1γ2 complex embedded in a 1-palmitoyl-2-oleoyl-phosphatidylcholine bilayer, using cross-linking information as validation. Our results show that although molecular dynamics simulations started with the G protein orientation in the β2-AR*·Gαsβ1γ2 complex crystal structure, the Gαi1β1γ2 protein reoriented itself within 300 ns. Two major changes occurred as follows. 1) The Gαi1 α5 helix tilt changed due to the outward movement of TMH5 in CB2 R*. 2) A 25° clockwise rotation of Gαi1β1γ2 underneath CB2 R* occurred, with rotation ceasing when Pro-139 (IC-2 loop) anchors in a hydrophobic pocket on Gαi1 (Val-34, Leu-194, Phe-196, Phe-336, Thr-340, Ile-343, and Ile-344). In this complex, all three experimentally identified cross-links can occur. These findings should be relevant for other class A G protein-coupled receptors that couple to Gi proteins.

The development of novel antibiotic drugs is one of the most pressing biomedical problems due to the increasing number of antibiotic-resistant pathogens. Antimicrobial peptides and lipopeptides are a promising category of candidates, but the molecular origins of their antimembrane activity is unclear. Here we explore a series of recently developed antimicrobial lipopeptides, using coarse-grained molecular-dynamics simulations and free energy methods to uncover the thermodynamics governing their binding to membranes. Specifically, we quantify C16-KGGK’s binding affinity to the two types of membrane by umbrella sampling. We also examined the origin of C16-KGGK’s selectivity for bacterial versus mammalian membranes by systematically varying the peptide sequence and salt concentration. Our data showed that the C16 hydrophobic tail is the main contributor to its affinity to lipid membrane, whereas the peptide portion is mainly responsible for its selectivity. Furthermore, the electrostatic interaction between the cationic peptide and anionic bacterial membrane plays a significant role in the selectivity.

G protein-coupled receptors (GPCRs) are a vital class of proteins that transduce biological signals across the cell membrane. However, their allosteric activation mechanism is not fully understood; crystal structures of active and inactive receptors have been reported, but the functional pathway between these two states remains elusive. Here, we employ structure-based (Gō-like) models to simulate activation of two GPCRs, rhodopsin and the β2 adrenergic receptor (β2AR). We used data-derived reaction coordinates that capture the activation mechanism for both proteins, showing that activation proceeds through quantitatively different paths in the two systems. Both reaction coordinates are determined from the dominant concerted motions in the simulations so the technique is broadly applicable. There were two surprising results. First, the main structural changes in the simulations were distributed throughout the transmembrane bundle, and not localized to the obvious areas of interest, such as the intracellular portion of helix 6. Second, the activation (and deactivation) paths were distinctly non-monotonic, populating states that were not simply interpolations between the inactive and active structures. These transitions also suggest a functional explanation for β2AR’s basal activity: it can proceed through a more broadly defined path during the observed transitions.

Several experimental structures of integral membrane proteins have demonstrated that tryptophan side chains are found preferentially within the membrane interface, but no clear physical explanation for this preference has been demonstrated. The present work contains a series of molecular dynamics simulations designed to provide insight into this phenomenon. Specifically, we performed conventional unbiased simulations, as well as biased molecular dynamics of indole and N-methylindole within palmitoyloleoyl phosphatidylcholine (POPC) bilayers. The unbiased calculations examined the behavior of the indoles in three environmentsthe membrane center, the headgroup region, and the water outside the membrane. The biased calculations were performed using a new dynamic windowing procedure, which efficiently moved the indoles through the bilayer. These calculations allowed us to explore the orientation, hydrogen bonding, and energetics of the indole, as a function of its location along the membrane normal. Furthermo...

A dipole lattice model for lipid membranes and their interactions with peptides is presented. It uses the Langevin dipole method to calculate electrostatic interactions in the heteroge- neous membrane environment. A series of test cases are presented, including spherical charges, dipoles, side chain analogs, and helical peptides. The model consistently produces qualitatively correct results. Proteins 2000;41:211-223. © 2000 Wiley-Liss, Inc.

The GPCR signaling cascade via the G protein pathway begins with an agonist binding to an inactive GPCR, causing conformational changes that activate the GPCR. Previously, we used molecular dynamics simulations to study the activation of the CB2 receptor (a Class A GPCR), by the endogenous ligand, 2-arachidonoylglycerol (2-AG) via the lipid bilayer (Hurst et al., 2010). In the next step of our study of the G-protein signaling cascade, we studied G-protein activation by an activated GPCR. In this step, GDP is released via the separation of the Gαi1 ras-like and helical domains of (Van Eps et al., 2011). To this end, we used our 2-AG activated CB2 model to produce an initial 2-AG/CB2/Gαi1β1γ2 complex based on the crystal structure of β2 adrenoreceptor in complex with Gαsβ1γ2 (Rasmussen et al., 2011). The complex was immersed in a fully hydrated 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer and NPT NAMD2 (Phillips et al., 2005) molecular dynamics (MD) simulations were initiated for four different trajectories of this system. Results from our longest MD simulation (3μs) suggest that the C-terminal α5 helix of Gαi prefers a different orientation in the CB2 activated receptor relative to the orientation seen in the empty state (GDP-GTP less) β2 adrenoreceptor/Gαsβ1γ2 complex. Initial hydrophobic interactions between P139 on CB2 intracellular loop 2 (IC-2 loop) and a hydrophobic pocket on Gαi consisting of residues V34 (N terminus); F336, T340, I343 and I344 (C terminus); L194 (β1-sheet); and F196 (β2-sheet) stabilize the receptor / Gαi1 protein interface during the first 0.5μs of the simulation. Later, between 1.4 - 1.6 μs, electrostatic interactions between R229 (CB2 IC-3 loop) and Q304/E308 (Gαi1 α4 helix) facilitate hydration of GDP. [Support: RO1 DA003934 and KO5 DA021358 (PHR)].

Protein function is a complicated interplay between structure and dynamics, which can be heavily influenced by environmental factors and conditions. This is particularly true in the case of membrane proteins, such as the visual receptor rhodopsin. It has been well documented that lipid headgroups, polyunsaturated tails, and the concentration of cholesterol in membranes all play a role in the function of rhodopsin. Recently, we used all-atom simulations to demonstrate that different lipid species have preferential interactions and possible binding sites on rhodopsin’s surface, consistent with experiment. However, the limited timescales of the simulations meant that the statistical uncertainty of these results was substantial. Accordingly, we present here 32 independent 1.6 µs coarse-grained simulations exploring lipids and cholesterols surrounding rhodopsin and opsin, in lipid bilayers mimicking those found naturally. Our results agree with those found experimentally and in previous simulations, but with far better statistical certainty. The results demonstrate the value of combining all-atom and coarse-grained models with experiment to provide a well-rounded view of lipid-protein interactions.

The emergence of antibiotic resistant pathogens is one of the major medical concerns of the 21st century, prompting renewed interest in the development of novel antimicrobial compounds. Here we use microsecond-scale all-atom molecular dynamics simulations to characterize the structure, dynamics, and membrane-binding mechanism of a synthetic antimicrobial lipopeptide, C16-KGGK. Our simulations suggest that these lipopeptides prefer to aggregate in solution and alter the intrinsic order of the lipid bilayer upon binding. From these results and previous coarse-grained simulations, we have developed a simple model for the binding and insertion process for these lipopeptides.

The global fold of human cannabinoid type 2 (CB2 ) receptor in the agonist-bound active state in lipid bilayers was investigated by solid-state (13)C- and (15)N magic-angle spinning (MAS) NMR, in combination with chemical-shift prediction from a structural model of the receptor obtained by microsecond-long molecular dynamics (MD) simulations. Uniformly (13)C- and (15)N-labeled CB2 receptor was expressed in milligram quantities by bacterial fermentation, purified, and functionally reconstituted into liposomes. (13)C MAS NMR spectra were recorded without sensitivity enhancement for direct comparison of Cα, Cβ, and C=O bands of superimposed resonances with predictions from protein structures generated by MD. The experimental NMR spectra matched the calculated spectra reasonably well indicating agreement of the global fold of the protein between experiment and simulations. In particular, the (13) C chemical shift distribution of Cα resonances was shown to be very sensitive to both the primary amino acid sequence and the secondary structure of CB2. Thus the shape of the Cα band can be used as an indicator of CB2 global fold. The prediction from MD simulations indicated that upon receptor activation a rather limited number of amino acid residues, mainly located in the extracellular Loop 2 and the second half of intracellular Loop 3, change their chemical shifts significantly (≥ 1.5 ppm for carbons and ≥ 5.0 ppm for nitrogens). Simulated two-dimensional (13) Cα(i)-(13)C=O(i) and (13)C=O(i)-(15)NH(i + 1) dipolar-interaction correlation spectra provide guidance for selective amino acid labeling and signal assignment schemes to study the molecular mechanism of activation of CB2 by solid-state MAS NMR.

The cannabinoid receptor subtype 2 (CB2) is a class A (rhodopsin-like) G-protein-coupled receptor (GPCR) that is highly expressed in the immune system (Galiegue et al., 1995) and signals via Gi proteins (Glass and Northup, 1999). The GPCR signaling cascade begins with agonist binding to an inactive receptor, causing conformational changes that activate the receptor. In previous work, we used molecular dynamics simulations to study the activation of the CB2 receptor, by the endogenous ligand, 2-arachidonoylglycerol (2-AG) via the lipid bilayer (Hurst et al., 2010). In work described here, we used our 2-AG activated CB2 model to produce an initial 2-AG/CB2/Gαi1β1γ2 assembly based on the crystal structure of β2 adrenoreceptor in complex with Gαsβ1γ2 (Rasmussen et al., 2011). Here, the G protein was located underneath CB2, but not inserted into CB2. The 2-AG/CB2/Gαi1β1γ2 assembly was then immersed in a fully hydrated 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer and NPT NAMD (Phillips et al., 2005) molecular dynamics simulations were initiated for four different makes of this system. Initial dynamics (50ns - 150ns) suggests that the C-terminus of Gαi (last 11 residues on α5 helix I344 to F354) inserts itself into the intracellular opening of the CB2 activated receptor at an altered angle compared to the β2 adrenoreceptor/Gαsβ1γ2 complex. The distance between intracellular ends of CB2 TMH3 and TMH6 increases relative to the 2-AG/CB2 complex, suggesting the first step in the formation of the 2-AG/CB2/Gαi1β1γ2 complex. [Support: RO1 DA003934 and KO5 DA021358 (PHR)]

LfB6 (RRWQWR-NH2) is a tryptophan- and arginine-rich cationic antimicrobial peptide with broad spectrum activity derived from bovine lactoferrin. Membrane binding occurs via electrostatic interactions between arginines and negative charges on the bacterial cell membrane and intercalation of the tryptophans at the membrane interface. N-terminal acylation (CH3(CH2)4CO-RRWQWR-NH2; C6-LfB6) can enhance the antimicrobial activity (Greathouse et al. (2008) J. Pept. Sci. 14:1103). Solid-state 2H and 31P NMR spectroscopy combined with all-atom and coarse-grained molecular dynamics (CG-MD) simulations have confirmed subtle differences between 1 mol% LfB6 and C6-LfB in bilayers composed of 3:1 POPE:POPG (anionic, bacterial-like) and POPC (zwitterionic, mammalian-like). MD simulations reveal that the arginines of C6-LfB6 make first contact with POPE:POPG; whereas the C6 tails are first to contact POPC. LfB6 shows no sequence preference. Additionally, C6-LfB6 inserts more deeply than LfB6 into both membranes. Tryptophan emission fluorescence spectra suggest the tryptophans in LfB6 and C6-LfB6 are more water exposed in neutral compared to anionic membranes, while CG-MD simulations reveal that LfB6 comes off the POPC membrane, exposing the tryptophans to water. Acylation, therefore, increases the “stickiness” of the peptide for lipid bilayers. Although both peptides at 1 mol% show significant membrane effects during short range simulations, C6-LfB6 has less influence on lipid order. We now compare experimental and molecular dynamics results for LfB6 and C6- LfB6 at 4 mol%. Solid-state 2H NMR spectra indicate that C6-LfB6 has a greater effect on the lipid acyl chain order at 4 mol% compared to 1 mol%; whereas the effects of LfB6 are similar at both concentrations. Molecular dynamics simulations will be presented for comparison.

HIV-1 reverse transcriptase is a critical drug target for HIV treatment, and understanding the exact mechanisms of its function and inhibition would significantly accelerate the development of new anti-HIV drugs. RT is a heterodimeric, multifunctional, multidomain protein with a 66 kDa subunit containing all of the catalytic active sites, and a 51 kDa subunit which is thought to provide structural stability to the larger subunit. Structural information on reverse transcriptase alone has proven to be insufficient to explain the mechanism of inhibition and drug resistance of non-nucleoside reverse transcriptase inhibitors. Elastic network modeling provides a technique to rapidly probe and compare protein dynamics. Combining elastic network modeling with hierarchical clusters of both structural and dynamic data reveals a wealth of novel information. Here we present an extensive survey of the dynamics of reverse transcriptase bound to a variety of ligands with a number of mutations, revealing a novel mechanism for drug resistance to non-nucleoside reverse transcriptase inhibitors, where hydrophobic core mutations subtly shift the position of the thumb subdomain, restoring active-state motion to multiple functionally significant regions of HIV-1 RT. This model arises out of a combination of structural and dynamic information, rather than exclusively from one or the other.

HIV-1 reverse transcriptase is a critical drug target for HIV treatment, and understanding the exact mechanisms of its function and inhibition would significantly accelerate the development of new anti-HIV drugs. It is well known that structure plays a critical role in protein function, but for reverse transcriptase, structural information has proven to be insufficient – despite enormous effort – to explain the mechanism of inhibition and drug resistance of non-nucleoside reverse transcriptase inhibitors. We hypothesize that the missing link is dynamics, information about the motions of the system. However, many of the techniques that give the best information about dynamics, such as solution NMR and molecular dynamics simulations, cannot be easily applied to a protein as large as reverse transcriptase. As an alternative, we combine elastic network modeling with simultaneous hierarchical clustering of structural and dynamic data. We present an extensive survey of the dynamics of reverse transcriptase bound to a variety of ligands and with a number of mutations, revealing a novel mechanism for drug resistance to non-nucleoside reverse transcriptase inhibitors. Hydrophobic core mutations restore active-state motion to multiple functionally significant regions of HIV-1 RT. This model arises out of a combination of structural and dynamic information, rather than exclusively from one or the other.

A series of synthetic antimicrobial lipopeptides (AMLPs) based around a common architecture of 4 amino acids (2 lysines), with a saturated fatty acid conjugated to the N-terminus, have been shown to have broad-spectrum antimicrobial activity and low hemolytic activity. Previous all-atom and coarse-grained molecular dynamics simulations from our group have shown that these molecules form micelles in solution and readily bind to model lipid bilayers. Here, we used microsecond-scale coarse-grained molecular dynamic simulations with the MARTINI forcefield to explore the thermodynamics governing the binding process, considering both isolated lipopeptides molecules and the micellar state. Using a combination of equilibrium umbrella sampling and non-equilibrium Jarzynski-style calculations, we estimate the binding free energy and explore the mechanism of entry. Our results provide biophysical insights into the mechanism of lipopeptides’ antimicrobial action.

The rhodopsin proteolipid system was modeled through all-atom molecular dynamics (MD) simulations on the microsecond timescale and compared to experimental 2H NMR data. The post-isomerization behavior of the covalently bound ligand, retinal, was tested for two independent models: (i) the counterion-switch (neutral binding pocket) and (ii) complex-counterion (negative binding pocket) [1]. Each model leads to distinct geometrical rearrangements, with direct implications for interpretation of recent rhodopsin crystal structures [2,3]. The distinctive feature of the counterion-switch simulation entails fluctuations of the C11=C12-C13=C14 dihedral angle. These motions are correlated with changes in the C5-, C9-, and C13-methyl group orientations, and result in a long-axis flip of the polyene chain within the binding pocket of rhodopsin. In contrast, the complex-counterion simulation produced retinal fluctuations at the C7=C8-C9=C10 dihedral. This facilitates stabilization of the methyl group orientations (≈60° with respect to the membrane normal), consistent with 2H NMR results for the Meta I state [4]. Recently, a putative structure of the fully-activated Meta II state revealed a long-axis flip of the retinylidene chain relative to its orientation in the dark state [2]. Our simulations show that electrostatic changes to the rhodopsin binding pocket lead to alternate pathways of retinal conformational release. Agreement between simulation and spectroscopy indicates that the retinylidene flip may only occur in the Meta I ↔ Meta II transition. These simulations provide an essential framework for interpreting molecular snapshots (crystal structures) of membrane protein activation, illuminating how small-scale changes in a GPCR binding pocket can affect large-scale membrane protein-lipid bilayer dynamics.[1] K. Martinez-Mayorga et al. (2006) JACS 128:16502.[2] H.-W. Choe et al. (2011) Nature 471:651.[3] J. Standfuss et al. (2011) Nature 471:656.[4] G.F.J. Salgado et al. (2006) JACS 128:11067.

The prevalence of antibiotic-resistant pathogens is a major medical concern, prompting increased interest in the development of novel antimicrobial compounds. One such set of naturally occurring compounds, known as antimicrobial peptides (AMPs), have broad-spectrum activity, but come with many limitations for clinical use. Recent work has resulted in a set of antimicrobial lipopeptides (AMLPs) with micromolar minimum inhibitory concentrations and excellent selectivity for bacterial membranes. To characterize a potent, synthetic lipopeptide, C16-KGGK, we used multi-microsecond coarse-grained simulations with the MARTINI forcefield, with a total simulation time of nearly 46μs. These simulations show rapid binding of C16-KGGK, which forms micelles in solution, to model bacterial lipid bilayers. Furthermore, upon binding to the surface of the bilayer, these lipopeptides alter the local lipid organization by recruiting negatively charged POPG lipids to the site of binding. It is likely that this drastic reorganization of the bilayer has major effects on bilayer dynamics and cellular processes that depend on specific bilayer compositions. By contrast, the simulations revealed no association between the lipopeptides and model mammalian bilayers. These simulations provide biophysical insights into lipopeptide selectivity and suggest a possible mechanism for antimicrobial action. This article is part of a Special Issue entitled: Membrane protein structure and function.

Elastic network models (ENMs) offer the capability to study collective motions of large systems without the overwhelming computational cost of all-atom molecular dynamics (MD). They accomplish this by replacing the standard force field terms with simple harmonic connections, typically between neighboring alpha carbons; the collective motions can then be computed using normal mode analysis. We investigate new methods for performing ENM calculations to improve their sensitivity and fidelity. First, we have explored using a higher resolution representation for the protein, supplementing the standard 1-node per residue approach by adding a second node to represent the side chain. Our aim is to utilize this increase in spatial resolution, enabling the ENMs to predict the effects of subtler changes, such as the binding of small ligands. Second, we have explored the efficacy of different choices for the spring constants in the harmonic representation. Following our previous work (Romo and Grossfield. Proteins, 2010), we have implemented a new functional form, which distinguishes bonded and nonbonded interactions. In both cases, we have validated our approaches via quantitative comparison to several microsecond-scale all-atom simulations of G protein-coupled receptors.

We investigated retinal conformational changes within the binding cavity of rhodopsin through all-atom molecular dynamics (MD) simulations on the microsecond time scale. Two distinct pathways involving post-isomerization release of retinal conformational strain were tested. The MD trajectories for the counterion-switch and complex-counterion mechanisms [1] show that retinal experiences completely distinct geometrical rearrangements, yielding differences in its orientation and conformation. In the counterion-switch simulation the dihedral angle C11=C12-C13=C14 fluctuates, and is correlated with changes in the C5-, C9-, and C13-methyl group orientations with respect to the membrane normal. For the complex-counterion simulation, changes in the dihedral angle C7=C8-C9=C10 allow the methyl group orientation to remain relatively unaffected (∼60° for the C5-, C9-, and C13-methyl groups), consistent with 2H NMR data [2,3]. The C5=C6-C7=C8 and C9=C10-C11=C12 torsion angles also experience significant fluctuations, but at different times in each counterion simulation and not in conjunction with retinal methyl re-orientation. In addition to retinal polyene chain analysis, the calculated 2H NMR spectra from each simulation clearly demonstrate that the retinal methyl orientations for the complex-counterion mechanism accurately reproduce experimental 2H NMR data for Meta I [2,3], whereas the counterion simulation does not. Rhodopsin simulations on the microsecond time scale are essential for modeling rhodopsin activation up to the Meta I photointermediate, because retinal polyene chain fluctuations occur as late as ∼1500 ns into each simulation. This work builds on our previous computational studies examining rhodopsin in the dark and Meta I states [1,4] and provides valuable insights into the rhodopsin activation process. [1] K. Martinez-Mayorga et al. (2006) JACS128, 16502-16503. [2] G.F.J. Salgado et al. (2006) JACS128, 11067-11071. [3] A.V. Struts et al. (2007) JMB372, 50-66. [4] P.-W. Lau et al. (2007) JMB372, 906-917.

The prevalence of antibiotic resistant pathogens is a major medical concern, prompting increased interest in the development of novel antimicrobial compounds. To characterize a potent, synthetic lipopeptide, C16-KGGK, microsecond time-scale all-atom simulations with the CHARMM forcefield are utilized. This lipopeptide targets the bacterial membrane, but the binding and reorganization processes are extremely slow. To increase simulation speed, supplemental coarse-grained simulations with the MARTINI forcefield are used. This combination provides insights into the structure, dynamics, binding and mechanism of antimicrobial action.View Large Image | View Hi-Res Image | Download PowerPoint Slide

In this issue, Neale et al. (1) present a calculation of the free energy to bind an antimicrobial peptide to a lipid bilayer using molecular dynamics simulations. This in itself is not unusual: many groups have used simulations to explore similar systems, and several have attempted to derive the binding thermodynamics. What is exceptional (and disturbing) about this article is the sheer computational effort required to get a good answer. Although Neale et al. (2) use a state-of-the-art Hamiltonian replica exchange technique, their results show that equilibration requires an astonishing 4 μs per simulation window. Worse yet, the results show that the error is not randomly distributed. Rather, the estimated free energy of binding becomes systematically more favorable as the runs are extended, suggesting that what we are seeing is an elongated relaxation process as opposed to simple improvements in statistical accuracy. These last two concepts are often conflated, but long relaxation times can cause quite different symptoms in a simulation from simple statistical error. This is best understood by considering the expected value of some property A〈y〉 computed from the simulation. If the main concern is simple statistical uncertainty, then we know two things: 1. As simulation time increases, 〈A〉 → Ao, where Ao is the true value (given the force field and simulation conditions); and 2. We can expect that if we run multiple trajectories, the values of the 〈A〉 value computed from the trajectory will be distributed randomly about Ao, with a variance that drops roughly as 1/Tsim, where Tsim is the length of the simulations. By contrast, systems with slow relaxations built in may not behave in this manner. For example, consider the system described by Fig. 1, which has two pairs of energy minima; each pair is separated by a small barrier in y, but a large barrier in x. If the quantity we are interested in is primarily a function of y (e.g. A(y)), the system will appear to make many transitions and 〈A〉 will appear to converge rapidly. However, if in building the system we consistently start on the left side of the x barrier, 〈A〉 will initially not converge to Ao. Rather, it will converge toward some different value A′o, representing averaging over the left-hand side of the conformation space. Moreover, the standard tools developed for examining a scalar time series, like autocorrelation analysis and block averaging (3), will fail to detect the problem, because the kinetics of y appear fast. Even more sophisticated global sampling assessment methods (4–6) may struggle to recognize the problem, particularly if no slow transitions occur at all; none of these methods can tell you what has not yet been seen. As a result, they are quite good at identifying mediocre to poor sampling, but less effective when the sampling is truly awful. Figure 1 Energy surface with fast and slow relaxing degrees of freedom. Although the kinetics in the y dimension appear fast, correct averaging of y will depend on sampling the x dimension as well, which contains a larger barrier that will cause slow sampling. ... Systems with slow relaxation will also behave differently as the simulation time Tsim is increased. Initially, the apparent uncertainty in A(y) will drop, but at longer times the variance will increase again, as the systems run long enough to occasionally cross the barrier in x; it is only when a large number of barrier crossings (and their reverse) have occurred that the 〈A〉 will converge to Ao. Ordinarily, the gold standard for quantifying error is to repeat the whole calculation using a different starting structure, but this too can fail if the construction procedure systematically produces one of the two states (e.g., you always start on the left side). Such circumstances are easily imaginable; for example, a crystal structure might capture one of two possible protein conformations. The case described by Neale et al. (1)—a peptide interacting with the membrane-water interface—represents another. It may sound like better system construction protocols will resolve these problems, and in principle they could. However, in practice, the fact is that the information needed to make optimal choices when building the system is generally not available, because that information is precisely what we hope to learn from the simulation. That molecular dynamics runs into serious challenges trying to obtain adequate sampling is hardly surprising—the biomolecular simulation field has a long history of undue optimism with respect to the timescales needed to get reliably interpretable results. What makes the work of Neale et al. (1) so impressive is that instead of running from the problem, they embraced it, and applied overwhelming computational effort to carefully characterize just how hard it is. It appears that there is no substitute for the investigators’ chemical intuition and sense of caution to identify in advance the likely timescales for structural transitions in the system. The publication of this work is a cautionary lesson to the rest of the simulation community about just how challenging this class of calculation is, and the kind of computational investment likely required to solve it.

Growing computing capacity and algorithmic advances have facilitated the study of increasingly large biomolecular systems at longer timescales. However, with these larger, more complex systems come questions about the quality of sampling and statistical convergence. What size systems can be sampled fully? If a system is not fully sampled, can certain “fast variables” be considered well-converged? How can one determine the statistical significance of observed results? The present review describes statistical tools and the underlying physical ideas necessary to address these questions. Basic definitions and ready-to-use analyses are provided, along with explicit recommendations. Such statistical analyses are of paramount importance in establishing the reliability of simulation data in any given study.