Thomas B. Woolf

MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations

MDAnalysis is an object‐oriented library for structural and temporal analysis of molecular dynamics (MD) simulation trajectories and individual protein structures. It is written in the Python language with some performance‐critical code in C. It uses the powerful NumPy package to expose trajectory data as fast and efficient NumPy arrays. It has been tested on systems of millions of particles. Many common file formats of simulation packages including CHARMM, Gromacs, Amber, and NAMD and the Protein Data Bank format can be read and written. Atoms can be selected with a syntax similar to CHARMMs powerful selection commands. MDAnalysis enables both novice and experienced programmers to rapidly write their own analytical tools and access data stored in trajectories in an easily accessible manner that facilitates interactive explorative analysis. MDAnalysis has been tested on and works for most Unix‐based platforms such as Linux and Mac OS X. It is freely available under the GNU General Public License from http://mdanalysis.googlecode.com.

Structure, energetics, and dynamics of lipid–protein interactions: A molecular dynamics study of the gramicidin A channel in a DMPC bilayer

The microscopic details of lipid–protein interactions are examined using molecular dynamics simulations of the gramicidin A channel embedded in a fully hydrated dimyristoyl phosphatidylcholine (DMPC) bilayer. A novel construction protocol was used to assemble the initial configurations of the membrane protein complex for the simulations. Three hundreds systems were constructed with different initial lipid placement and conformations. Seven systems were simulated with molecular dynamics. One system was simulated for a total of 600 psec, four were simulated for 300 psec, and two for 100 psec. Analysis of the resulting trajectories shows that the bulk solvent–membrane interface region is much broader than traditionally pictured in simplified continuum theories: its width is almost 15 Å. In addition, lipid–protein interactions are far more varied, both structurally and energetically, than is usually assumed: the total interaction energy between the gramicidin A and the individual lipids varies from 0 to −50 kcal/mol. The deuterium quadrupolar splittings of the lipid acyl chains calculated from the trajectories are in good agreement with experimental data. The lipid chains in direct contact with the GA are ordered but the effect is not uniform due to the irregular surface of the protein. Energy decompositions shows that the most energetically favorable interactions between lipid and protein involve nearly equal contributions from van der Waals and electrostatic interactions. The tryptophans, located near the bulk–membrane interface, appear to be particularly important in mediating both hydrogen bonding interactions with the lipid glycerol backbone and water and also in forming favorable van der Waals contacts with the hydrocarbon chains. In contrast, the interactions of the leucine residues with the lipids, also located near the interface, are dominated by van der Waals interactions with the hydrocarbon lipid chains.

Atomistic simulations of biologically realistic transmembrane potential gradients.

We present all-atom molecular dynamics simulations of biologically realistic transmembrane potential gradients across a DMPC bilayer. These simulations are the first to model this gradient in all-atom detail, with the field generated solely by explicit ion dynamics. Unlike traditional bilayer simulations that have one bilayer per unit cell, we simulate a 170 mV potential gradient by using a unit cell consisting of three salt-water baths separated by two bilayers, with full three-dimensional periodicity. The study shows that current computational resources are powerful enough to generate a truly electrified interface, as we show the predicted effect of the field on the overall charge distribution. Additionally, starting from Poissons equation, we show a new derivation of the double integral equation for calculating the potential profile in systems with this type of periodicity.

Improving the efficiency and reliability of free energy perturbation calculations using overlap sampling methods

A challenge in free energy calculation for complex molecular systems by computer simulation is to obtain a reliable estimate within feasible computational time. In this study, we suggest an answer to this challenge by exploring a simple method, overlap sampling (OS), for producing reliable free‐energy results in an efficient way. The formalism of the OS method is based on ensuring sampling of important overlapping phase space during perturbation calculations. This technique samples both forward and reverse free energy perturbation (FEP) to improve the free‐energy calculation. It considers the asymmetry of the FEP calculation and features an ability to optimize both the precision and the accuracy of the measurement without affecting the simulation process itself. The OS method is tested at two optimization levels: no optimization (simple OS), and full optimization (equivalent to Bennetts method), and compared to conventional FEP techniques, including the widely used direct FEP averaging method, on three alchemical mutation systems: (a) an anion transformation in water solution, (b) mutation between methanol and ethane, and (c) alchemical change of an adenosine molecule. It is consistently shown that the reliability of free‐energy estimates can be greatly improved using the OS techniques at both optimization levels, while the performance of Bennetts method is particularly striking. In addition, the efficiency of a calculation can be significantly improved because the method is able to (a) converge to the right answer quickly, and (b) work for large perturbations. The basic two‐stage OS method can be extended to admit additional stages, if needed. We suggest that the OS method can be used as a general perturbation technique for computing free energy differences in molecular simulations.

Reduced amino acid alphabet is sufficient to accurately recognize intrinsically disordered protein

Intrinsically disordered proteins are an important class of proteins with unique functions and properties. Here, we have applied a support vector machine (SVM) trained on naturally occurring disordered and ordered proteins to examine the contribution of various parameters (vectors) to recognizing proteins that contain disordered regions. We find that a SVM that incorporates only amino acid composition has a recognition accuracy of 87 ± 2%. This result suggests that composition alone is sufficient to accurately recognize disorder. Interestingly, SVMs using reduced sets of amino acids based on chemical similarity preserve high recognition accuracy. A set as small as four retains an accuracy of 84 ± 2%; this suggests that general physicochemical properties rather than specific amino acids are important factors contributing to protein disorder.

Manganese Selectivity of Pmr1, the Yeast Secretory Pathway Ion Pump, Is Defined by Residue Gln783 in Transmembrane Segment 6 RESIDUE Asp778 IS ESSENTIAL FOR CATION TRANSPORT

We have solubilized and purified the histidine-tagged yeast secretory pathway/Golgi ion pump Pmr1 to near homogeneity in one step, using nickel affinity chromatography. The purified pump demonstrates both Ca2+- and Mn2+-dependent ATP hydrolysis and phosphoenzyme intermediate formation in forward (ATP) and reverse (Pi) directions. This preparation has allowed us to examine, in detail, the properties of mutations D778A and Q783A in transmembrane segment M6 of Pmr1. In phenotypic screens of Ca2+ chelator and Mn2+ toxicity reported separately (Wei, Y., Chen, J., Rosas, G., Tompkins, D.A., Holt, P.A., and Rao, R. (2000) J. Biol. Chem. 275, XXXX–XXXX), D778A was a loss-of-function mutant apparently defective for transport of both Ca2+ and Mn2+, whereas mutant Q783A displayed a differential sensitivity consistent with the selective loss of Mn2+transport. We show that mutant D778A is devoid of cation-dependent ATP hydrolytic activity and phosphoenzyme formation from ATP. However, reverse phosphorylation from Pi is preserved but is insensitive to inhibition by Ca2+ or Mn2+ ions, which is evidence for a specific inability to bind cations in this mutant. We also show that Ca2+ can activate ATP hydrolysis in the purified Q783A mutant, with a half-maximal concentration of 0.06 μm, essentially identical to that of wild type (0.07 μm). Mn2+ activation of ATP hydrolysis was half-maximal at 0.02 μm in wild type, establishing a normal selectivity profile of Mn2+ > Ca2+. Strikingly, Mn2+-ATPase in the Q783A mutant was nearly abolished, even at concentrations of up to 10 μm. These results were confirmed in assays of phosphoenzyme intermediates. Molecular modeling of the packing between helices M4 and M6 suggests that residue Gln783 in M6 may form a critical hydrophobic interaction with Val335 in M4, such that the Ala substitution modifies the packing or tilt of the helices and thus the ion pore. The data emphasize the critical role of transmembrane segment M6 in defining the cation binding pocket of P-type ATPases.

Folding is not required for bilayer insertion: Replica exchange simulations of an α‐helical peptide with an explicit lipid bilayer

We implement the replica exchange molecular dynamics algorithm to study the interactions of a model peptide (WALP‐16) with an explicitly represented DPPC membrane bilayer. We observe the spontaneous, unbiased insertion of WALP‐16 into the DPPC bilayer and its folding into an α‐helix with a transbilayer orientation. The free energy surface suggests that the insertion of the peptide into the DPPC bilayer precedes secondary structure formation. Although the peptide has some propensity to form a partially helical structure in the interfacial region of the DPPC/water system, this state is not a productive intermediate but rather an off‐pathway trap for WALP‐16 insertion. Equilibrium simulations show that the observed insertion/folding pathway mirrors the potential of mean force (PMF). Calculation of the enthalpic and entropic contributions to this PMF show that the surface bound conformation of WALP‐16 is significantly lower in energy than other conformations, and that the insertion of WALP‐16 into the bilayer without regular secondary structure is enthalpically unfavorable by 5–10 kcal/mol/residue. The observed insertion/folding pathway disagrees with the dominant conceptual model, 1–3 which is that a surface‐bound helix is an obligatory intermediate for the insertion of α‐helical peptides into lipid bilayers. In our simulations, the observed insertion/folding pathway is favored because of a large (>100 kcal/mol) increase in system entropy that occurs when the unstructured WALP‐16 peptide enters the lipid bilayer interior. The insertion/folding pathway that is lowest in free energy depends sensitively on the near cancellation of large enthalpic and entropic terms. This suggests the possibility that intrinsic membrane peptides may have a diversity of insertion/folding behaviors depending on the exact system of peptide and lipid under consideration. Proteins 2005.

Dynamic reaction paths and rates through importance-sampled stochastic dynamics

We extend a previously developed method, based on Wagner’s stochastic formulation of importance sampling, to the calculation of reaction rates and to a simple quantitative description of finite-temperature, average dynamic paths. Only the initial and final states are required as input—no information on transition state(s) is necessary. We demonstrate the method for a single particle moving on the two-dimensional Muller–Brown potential surface. Beyond computing the forward and reverse rates for this surface, we determine the average path, which exhibits “saddle point avoidance.” The method may be generalized to arbitrary numbers of degrees of freedom and to arbitrary types of stochastic dynamics.

The binding site of sodium in the gramicidin A channel: comparison of molecular dynamics with solid-state NMR data

The location of the main binding site for sodium in the gramicidin A (GA) channel was investigated with molecular dynamics simulations, using an atomic model of the channel embedded in a fully hydrated dimyristoyl phosphatidycholine (DMPC) bilayer. Twenty-four separate simulations in which a sodium was restrained at different locations along the channel axis were generated. The results are compared with carbonyl 13C chemical shift anisotropy solid-state NMR experimental data previously obtained with oriented GA:DMPC samples. Predictions are made for other solid-state NMR properties that could be observed experimentally. The combined information from experiment and simulation strongly suggests that the main binding sites for sodium are near the channels mouth, approximately 9.2 A from the center of the dimer channel. The 13C chemical shift anisotropy of Leu10 is the most affected by the presence of a sodium ion in the binding site. In the binding site, the sodium ion is lying off-axis, making contact with two carbonyl oxygens and two single-file water molecules. The main channel ligand is provided by the carbonyl group of the Leu10-Trp11 peptide linkage, which exhibits the largest deviation from the ion-free channel structure. Transient contacts with the carbonyl group of Val8 and Trp15 are also present. The influence of the tryptophan side chains on the channel conductance is examined based on the current information about the binding site.

Structure and dynamics of an amphiphilic peptide in a lipid bilayer: a molecular dynamics study

A molecular dynamics simulation of a simple model membrane system composed of a single amphiphilic helical peptide (ace-K2GL16K2A-amide) in a fully hydrated 1,2-dimyristoyl-sn-glycero-3-phosphocholine bilayer was performed for a total of 1060 ps. The secondary structure of the peptide and its stability were described in terms of average dihedral angles, phi and psi, and the C alpha torsion angles formed by backbone atoms; by the average translation per residue along the helix axis; and by the intramolecular peptide hydrogen bonds. The results indicated that residues 6 through 15 remain in a stable right-handed alpha-helical conformation, whereas both termini exhibit substantial fluctuations. A change in the backbone dihedral angles for residues 16 and 17 is accompanied by the loss of two intramolecular hydrogen bonds, leading to a local but long-lived disruption of the helix. The dynamics of the peptide was characterized in terms of local and global helix motions. The local motions of the N-H bond angles were described in terms of the autocorrelation functions of P2[cos thetaNH(t, t + tau)] and reflected the different degrees of local peptide order as well as a variation in time scale for local motions. The chi1 and chi2 dihedral angles of the leucine side chains underwent frequent transitions between potential minima. No connection between the side-chain positions and their mobility was observed, however. In contrast, the lysine side chains displayed little mobility during the simulation. The global peptide motions were characterized by the tilting and bending motions of the helix. Although the peptide was initially aligned parallel to the bilayer normal, during the simulation it was observed to tilt away from the normal, reaching an angle of approximately 25 degrees by the end of the simulation. In addition, a slight bend of the helix was detected. Finally, the solvation of the peptide backbone and side-chain atoms was also investigated.