Tod D. Romo

Validating and improving elastic network models with molecular dynamics simulations.

Elastic network models (ENMs) are a class of simple models intended to represent the collective motions of proteins. In contrast to all‐atom molecular dynamics simulations, the low computational investment required to use an ENM makes them ideal for speculative hypothesis‐testing situations. Historically, ENMs have been validated via comparison to crystallographic B‐factors, but this comparison is relatively low‐resolution and only tests the predictions of relative flexibility. In this work, we systematically validate and optimize a number of ENM‐type models by quantitatively comparing their predictions to microsecond‐scale all‐atom simulations of three different G protein coupled receptors. We show that, despite their apparent simplicity, well‐optimized ENMs perform remarkably well, reproducing the protein fluctuations with an accuracy comparable to what one would expect from all‐atom simulations run for several hundred nanoseconds. Proteins 2010.

Characterization of a potent antimicrobial lipopeptide via coarse-grained molecular dynamics

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.

Unknown Unknowns: the Challenge of Systematic and Statistical Error in Molecular Dynamics Simulations

In this issue, Neale etxa0al. (1) present axa0calculation of the free energy to bindxa0an 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 isxa0the sheer computational effort required to get a good answer. Although Neale etxa0al. (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, thexa0results 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. n n nThese 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. Ifxa0the main concern is simple statistical uncertainty, then we know two things: n n n1. n nAs simulation time increases, 〈A〉 → Ao, where Ao is the true value (given the force field and simulation conditions); and n n n2. n nWe 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. n n n n nBy contrast, systems with slow relaxations built in may not behave inxa0this manner. For example, consider the system described by Fig.xa01, 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 toxa0Ao. 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. n n n nFigure 1 n nEnergy 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. ... n n n nSystems 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. n nOrdinarily, 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 etxa0al. (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. n nThat molecular dynamics runs into serious challenges trying to obtain adequate sampling is hardly surprising—the biomolecular simulation fieldxa0has a long history of undue optimism with respect to the timescales needed to get reliably interpretable results. What makes the work of Neale etxa0al. (1) so impressive is that instead ofxa0running 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 thexa0likely 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 ofxa0calculation is, and the kind of computational investment likely required to solve it.

Membrane binding of an acyl-lactoferricin B antimicrobial peptide from solid-state NMR experiments and molecular dynamics simulations.

One approach to the growing health problem of antibiotic resistant bacteria is the development of antimicrobial peptides (AMPs) as alternative treatments. The mechanism by which these AMPs selectively attack the bacterial membrane is not well understood, but is believed to depend on differences in membrane lipid composition. N-acylation of the small amidated hexapeptide, RRWQWR-NH(2) (LfB6), derived from the 25 amino acid bovine lactoferricin (LfB25) can be an effective means to improve its antimicrobial properties. Here, we investigate the interactions of C6-LfB6, N-acylated with a 6 carbon fatty acid, with model lipid bilayers with two distinct compositions: 3:1 POPE:POPG (negatively charged) and POPC (zwitterionic). Results from solid-state (2)H and (31)P NMR experiments are compared with those from an ensemble of all-atom molecular dynamic simulations running in aggregate more than 8.6ms. (2)H NMR spectra reveal no change in the lipid acyl chain order when C6-LfB6 is bound to the negatively charged membrane and only a slight decrease in order when it is bound to the zwitterionic membrane. (31)P NMR spectra show no significant perturbation of the phosphate head groups of either lipid system in the presence of C6-LfB6. Molecular dynamic simulations show that for the negatively charged membrane, the peptides arginines drive the initial association with the membrane, followed by attachment of the tryptophans at the membrane-water interface, and finally by the insertion of the C6 tails deep into the bilayer. In contrast, the C6 tail leads the association with the zwitterionic membrane, with the tryptophans and arginines associating with the membrane-water interface in roughly the same amount of time. We find similar patterns in the order parameters from our simulations. Moreover, we find in the simulations that the C6 tail can insert 1-2Å more deeply into the zwitterionic membrane and can exist in a wider range of angles than in the negatively charged membrane. We propose this is due to the larger area per lipid in the zwitterionic membrane, which provides more space for the C6 to insert and assume different orientations.

Lightweight Object Oriented Structure Analysis: Tools for Building Tools to Analyze Molecular Dynamics Simulations

LOOS (Lightweight Object Oriented Structure‐analysis) is a C++ library designed to facilitate making novel tools for analyzing molecular dynamics simulations by abstracting out the repetitive tasks, allowing developers to focus on the scientifically relevant part of the problem. LOOS supports input using the native file formats of most common biomolecular simulation packages, including CHARMM, NAMD, Amber, Tinker, and Gromacs. A dynamic atom selection language based on the C expression syntax is included and is easily accessible to the tool‐writer. In addition, LOOS is bundled with over 140 prebuilt tools, including suites of tools for analyzing simulation convergence, three‐dimensional histograms, and elastic network models. Through modern C++ design, LOOS is both simple to develop with (requiring knowledge of only four core classes and a few utility functions) and is easily extensible. A python interface to the core classes is also provided, further facilitating tool development.

Retinal Conformation Changes Rhodopsin’s Dynamic Ensemble

G protein-coupled receptors are vital membrane proteins that allosterically transduce biomolecular signals across the cell membrane. However, the process by which ligand binding induces protein conformation changes is not well understood biophysically. Rhodopsin, the mammalian dim-light receptor, is a unique test case for understanding these processes because of its switch-like activity; the ligand, retinal, is bound throughout the activation cycle, switching from inverse agonist to agonist after absorbing a photon. By contrast, the ligand-free opsin is outside the activation cycle and may behave differently. We find that retinal influences rhodopsin dynamics using an ensemble of all-atom molecular dynamics simulations that in aggregate contain 100 μs of sampling. Active retinal destabilizes the inactive state of the receptor, whereas the active ensemble was more structurally homogenous. By contrast, simulations of an active-like receptor without retinal present were much more heterogeneous than those containing retinal. These results suggest allosteric processes are more complicated than a ligand inducing protein conformational changes or simply capturing a shifted ensemble as outlined in classic models of allostery.

Structure-Based Simulations Reveal Concerted Dynamics of GPCR Activation

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 use 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 nonmonotonic, populating states that were not simply interpolations between the inactive and active structures. These transitions also suggest a functional explanation for β2ARs basal activity: it can proceed through a more broadly defined path during the observed transitions. Proteins 2014; 82:2538–2551.

How Fast is Your Camera? Timescales for Molecular Motion and their Role in Restraining Molecular Dynamics

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. n nIn this issue of Biophysical Journal, De Simone etxa0al. (1) describe their efforts to build structural models fromxa0this 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 toxa0drive 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). n nAlthough other researchers (particularly Im etxa0al. (2), De Simone etxa0al. (3), and Richter etxa0al. (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-statexa02H 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. n nThere 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. n nThe 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. n nHowever, 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). n nBy contrast, the time-resolution of a MD simulation is specified by the time step (generally 2xa0fs 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. n nThere are three distinct cases that determine what sort of effects this kind of averaging can cause, shown pictorially in Fig.xa01. n n n nFigure 1 n nTimescales 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. n n n nIn 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 etxa0al. (1) will likely be successful. n nThe 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, anxa0ensemble-based approach, such as that proposed here by De Simone etxa0al. (1) and elsewhere by Im etxa0al. (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. n nHowever, 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.xa01), 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 ofxa0structures from the molecular dynamics; one could imagine that the formalism would resemble self-guided molecular dynamics (14). n nDespite this reservation, it is clear that the work of De Simone etxa0al. (1) is of significant value. MD simulations and solid-state NMR are two of thexa0most powerful techniques for exploring the structure and dynamics of membrane proteins, so methods to couple them have the potential for great synergy.