Transmembrane helix assembly by window exchange umbrella sampling

GENESIS 2.1: High-Performance Molecular Dynamics Software for Enhanced Sampling and Free-Energy Calculations for Atomistic, Coarse-Grained, and Quantum Mechanics/Molecular Mechanics Models

GENeralized-Ensemble SImulation System (GENESIS) is a molecular dynamics (MD) software developed to simulate the conformational dynamics of a single biomolecule, as well as molecular interactions in large biomolecular assemblies and between multiple biomolecules in cellular environments. To achieve the latter purpose, the earlier versions of GENESIS emphasized high performance in atomistic MD simulations on massively parallel supercomputers, with or without graphics processing units (GPUs). Here, we implemented multiscale MD simulations that include atomistic, coarse-grained, and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations. They demonstrate high performance and are integrated with enhanced conformational sampling algorithms and free-energy calculations without using external programs except for the QM programs. In this article, we review new functions, molecular models, and other essential features in GENESIS version 2.1 and discuss ongoing developments for future releases.

Advanced computational approaches to understand protein aggregation

Protein aggregation is a widespread phenomenon implicated in debilitating diseases like Alzheimer's, Parkinson's, and cataracts, presenting complex hurdles for the field of molecular biology. In this review, we explore the evolving realm of computational methods and bioinformatics tools that have revolutionized our comprehension of protein aggregation. Beginning with a discussion of the multifaceted challenges associated with understanding this process and emphasizing the critical need for precise predictive tools, we highlight how computational techniques have become indispensable for understanding protein aggregation. We focus on molecular simulations, notably molecular dynamics (MD) simulations, spanning from atomistic to coarse-grained levels, which have emerged as pivotal tools in unraveling the complex dynamics governing protein aggregation in diseases such as cataracts, Alzheimer's, and Parkinson's. MD simulations provide microscopic insights into protein interactions and the subtleties of aggregation pathways, with advanced techniques like replica exchange molecular dynamics, Metadynamics (MetaD), and umbrella sampling enhancing our understanding by probing intricate energy landscapes and transition states. We delve into specific applications of MD simulations, elucidating the chaperone mechanism underlying cataract formation using Markov state modeling and the intricate pathways and interactions driving the toxic aggregate formation in Alzheimer's and Parkinson's disease. Transitioning we highlight how computational techniques, including bioinformatics, sequence analysis, structural data, machine learning algorithms, and artificial intelligence have become indispensable for predicting protein aggregation propensity and locating aggregation-prone regions within protein sequences. Throughout our exploration, we underscore the symbiotic relationship between computational approaches and empirical data, which has paved the way for potential therapeutic strategies against protein aggregation-related diseases. In conclusion, this review offers a comprehensive overview of advanced computational methodologies and bioinformatics tools that have catalyzed breakthroughs in unraveling the molecular basis of protein aggregation, with significant implications for clinical interventions, standing at the intersection of computational biology and experimental research.

The α C helix is a central regulator of PKR activation

Protein kinase R (PKR) functions in the eukaryotic innate immune system as a first-line defense against viral infections. PKR binds viral dsRNA, leading to autophosphorylation and activation. In its active state, PKR can phosphorylate its primary substrate, eIF2 α , which blocks initiation of translation in the infected cell. It has been established that PKR activation occurs when the kinase domain dimerizes in a back-to-back configuration. However, the mechanism by which dimerization leads to enzymatic activation is not fully understood. Here, we investigate the structural mechanistic basis and energy landscape for PKR activation, with a focus on the α C helix - a kinase activation and signal integration hub - using all-atom equilibrium and enhanced sampling molecular dynamics simulations. By employing window-exchange umbrella sampling, we compute free energy profiles of activation which show that back-to-back dimerization stabilizes a catalytically competent conformation of PKR. Key hydrophobic residues in the homodimer interface contribute to stabilization of the α C helix in an active conformation and the position of its glutamate residue. Using linear mutual information analysis, we analyze allosteric communication connecting the protomers' N-lobes and the α C helix dimer interface with the α C helix.

Computational Insights into the Adsorption of Ligands on Gold Nanosurfaces

We study the adsorption process of model peptides, nucleobases, and selected standard ligands on gold through the development of a computational protocol based on fully atomistic classical molecular dynamics (MD) simulations combined with umbrella sampling techniques. The specific features of the interface components, namely, the molecule, the metallic substrate, and the solvent, are taken into account through different combinations of force fields (FFs), which are found to strongly affect the results, especially changing absolute and relative adsorption free energies and trends. Overall, noncovalent interactions drive the process along the adsorption pathways. Our findings also show that a suitable choice of the FF combinations can shed light on the affinity, position, orientation, and dynamic fluctuations of the target molecule with respect to the surface. The proposed protocol may help the understanding of the adsorption process at the microscopic level and may drive the in-silico design of biosensors for detection purposes.

Unsaturated Lipids Facilitate Partitioning of Transmembrane Peptides into the Liquid Ordered Phase

The affinity of single-pass transmembrane (TM) proteins for ordered membrane phases has been reported to depend on their lipidation, TM length, and lipid accessible surface area. In this work, the raft affinities of the TM domain of the linker for activation of T cells and its depalmitoylated variant are assessed using free energy simulations in a binary bilayer system composed of two laterally patched bilayers of ternary liquid ordered (Lo) and liquid disordered (Ld) phases. These phases are modeled by distinct compositions of distearoylphosphatidylcholine, palmitoyloleoylphosphatidylcholine (POPC), and cholesterol, and the simulations were carried out for 4.5 μs/window. Both peptides are shown to preferentially partition into the Ld phase in agreement with model membrane experiments and previous simulations on ternary lipid mixtures but not with measurements on giant plasma membrane vesicles where the Lo is slightly preferred. However, the 500 ns average relaxation time of lipid rearrangement around the peptide precluded a quantitative analysis of free energy differences arising from peptide palmitoylation and two different lipid compositions. When in the Lo phase, peptides reside in regions rich in POPC and interact preferentially with its unsaturated tail. Hence, the detailed substructure of the Lo phase is an important modulator of peptide partitioning, in addition to the inherent properties of the peptide.

CHARMM-GUI Enhanced Sampler for various collective variables and enhanced sampling methods

Enhanced sampling methodologies modifying underlying Hamiltonians can be used for the systems with a rugged potential energy surface that makes it hard to observe convergence using conventional unbiased molecular dynamics (MD) simulations. We present CHARMM-GUI Enhanced Sampler, a web-based tool to prepare various enhanced sampling simulations inputs with user-selected collective variables (CVs). Enhanced Sampler provides inputs for the following nine methods: accelerated MD, Gaussian accelerated MD, conformational flooding, metadynamics, adaptive biasing force, steered MD, temperature replica exchange MD, replica exchange solute tempering 2, and replica exchange umbrella sampling for the method-implemented MD packages including AMBER, CHARMM, GENESIS, GROMACS, NAMD, and OpenMM. Users only need to select a group of atoms via intuitive web-implementation in order to define commonly used nine CVs of interest: center of mass based distance, angle, dihedral, root-mean-square-distance, radius of gyration, distance projected on axis, two types of angles projected on axis, and coordination numbers. The enhanced sampling methods are tested with several biological systems to illustrate their efficiency over conventional MD. Enhanced Sampler with carefully optimized system-dependent parameters will help users to get meaningful results from their enhanced sampling simulations.

Molecular simulation of lignin-related aromatic compound permeation through gram-negative bacterial outer membranes

Lignin, an abundant aromatic heteropolymer in secondary plant cell walls, is the single largest source of renewable aromatics in the biosphere. Leveraging this resource for renewable bioproducts through targeted microbial action depends on lignin fragment uptake by microbial hosts and subsequent enzymatic action to obtain the desired product. Recent computational work has emphasized that bacterial inner membranes are permeable to many aromatic compounds expected from lignin depolymerization processes. In this study, we expand on these findings through simulations for 42 lignin-related compounds across a gram-negative bacterial outer membrane model. Unbiased simulation trajectories indicate that spontaneous crossing for the full outer membrane is relatively rare at molecular simulation timescales, primarily due to preferential membrane partitioning and slow diffusion within the lipopolysaccharide layer within the outer membrane. Membrane partitioning and permeability coefficients were determined through replica exchange umbrella sampling simulations to overcome sampling limitations. We find that the glycosylated lipopolysaccharides found in the outer membrane increase the permeation barrier to many lignin-related compounds, particularly the most hydrophobic compounds. However, the effect is relatively modest; at industrially relevant concentrations, uncharged lignin-related compounds will readily diffuse across the outer membrane without the need for specific porins. Together, our results provide insight into the permeability of the bacterial outer membrane for assessing lignin fragment uptake and the future production of renewable bioproducts.

Iterative Landmark-Based Umbrella Sampling (ILBUS) Protocol for Sampling of Conformational Space of Biomolecules

Computer simulations of biomolecules such as molecular dynamics often suffer from insufficient sampling. Due to limited computational resources, insufficient sampling prevents obtaining proper equilibrium distributions of observed properties. To deal with this problem, we proposed a simulation protocol for efficient resampling of collected off-equilibrium trajectories. These trajectories are utilized for the initial mapping of the conformational space, which is later properly resampled by the introduced Iterative Landmark-Based Umbrella Sampling (ILBUS) method. Reconstruction of static equilibrium properties is achieved by the multistate Bennett acceptance ratio (MBAR) method, which enables efficient use of simulated data. The ILBUS protocol is geometry-based and does not demand any additional collective variable or a dimensional-reduction technique. The only requirement is a set of suitably spaced reference conformations, which serve as landmarks in the mapped conformational space. Additionally, the ILBUS protocol encompasses an iterative process that optimizes the force constant used in the umbrella sampling simulation. Such tuning is an inherent feature of the protocol and does not need to be performed by the user in advance. Furthermore, even the simulations with suboptimal force constants can be used in estimates by MBAR. We demonstrate the feasibility and the performance of this approach in the study of the conformational landscape of the alanine dipeptide, met-enkephalin, and adenylate kinase.

Developing initial conditions for simulations of asymmetric membranes: a practical recommendation

It has been proposed that the surface tension difference between leaflets (or differential stress) in asymmetric bilayers is generally nonvanishing. This implies that there is no unique approach to generate initial conditions for simulations of asymmetric bilayers in the absence of experimentally derived constraints. Current generation methods include individual area per lipid (APL) based, leaflet surface area (SA) matching, and zero leaflet tension based (0-DS). This work adds a bilayer-based approach that aims for achieving partial chemical equilibrium by interleaflet switching of selected lipids via P2₁ periodic boundary conditions. Based on a recently proposed theoretical framework, we obtained expressions for tensions in asymmetric bilayers from both the bending and area strains. We also developed a quantitative measure for the energetic penalty from the differential stress. The impacts of APL-, SA-, and 0-DS-based approaches on mechanical properties are assessed for two different asymmetric bilayers. The lateral pressure profile and its moments differ significantly for each method, whereas the area compressibility modulus is relatively insensitive. Application of P2₁ periodic boundary conditions (APL/P2₁, SA/P2₁, and 0-DS/P2₁) results in better agreement in mechanical properties between asymmetric bilayers generated by APL-, SA-, and 0-DS-based approaches, in which changes are the smallest for bilayers from the SA-based method. The estimated differential stress from the theory shows good agreement with that from the simulations. These simulation results and the good agreement between the predicted and observed differential stress further support the theoretical framework in which bilayer mechanical properties are outcomes of the interplay between intrinsic bending and asymmetric lipid packing. Based on the simulation results and theoretical predictions, the SA/P2₁-based, or at least the SA-based (when the differential stress is small), approach is recommended as a practical method for developing initial conditions for asymmetric bilayer simulations.

Vectorial insertion of a β-helical peptide into membrane: a theoretical study on polytheonamide B

Spontaneous unidirectional, or vectorial, insertion of transmembrane peptides is a fundamental biophysical process for toxin and viral actions. Polytheonamide B (pTB) is a potent cytotoxic peptide with a β^6.3-helical structure. Previous experimental studies revealed that the pTB inserts into the membrane in a vectorial fashion and forms a channel with its single molecular length long enough to span the membrane. Also, molecular dynamics simulation studies demonstrated that the pTB is prefolded in aqueous solution. These are unique features of pTB because most of the peptide toxins form channels through oligomerization of transmembrane helices. Here, we performed all-atom molecular dynamics simulations to examine the dynamic mechanism of the vectorial insertion of pTB, providing underlying elementary processes of the membrane insertion of a prefolded single transmembrane peptide. We find that the insertion of pTB proceeds with only the local lateral compression of the membrane in three successive phases: "landing," "penetration," and "equilibration" phases. The free energy calculations using the replica-exchange umbrella sampling simulations present an energy cost of 4.3 kcal/mol at the membrane surface for the membrane insertion of pTB from bulk water. The trajectories of membrane insertion revealed that the insertion process can occur in two possible pathways, namely "trapped" and "untrapped" insertions; in some cases, pTB is trapped in the upper leaflet during the penetration phase. Our simulations demonstrated the importance of membrane anchoring by the hydrophobic N-terminal blocking group in the landing phase, leading to subsequent vectorial insertion.

Conformational Fluctuations in GTP-Bound K-Ras: A Metadynamics Perspective with Harmonic Linear Discriminant Analysis

Biomacromolecules often undergo significant conformational rearrangements during function. In proteins, these motions typically consist in nontrivial, concerted rearrangement of multiple flexible regions. Mechanistic, thermodynamics, and kinetic predictions can be obtained via molecular dynamics simulations, provided that the simulation time is at least comparable to the relevant time scale of the process of interest. Because of the substantial computational cost, however, plain MD simulations often have difficulty in obtaining sufficient statistics for converged estimates, requiring the use of more-advanced techniques. Central in many enhanced sampling methods is the definition of a small set of relevant degrees of freedom (collective variables) that are able to describe the transitions between different metastable states of the system. The harmonic linear discriminant analysis (HLDA) has been shown to be useful for constructing low-dimensional collective variables in various complex systems. Here, we apply HLDA to study the free-energy landscape of a monomeric protein around its native state. More precisely, we study the K-Ras protein bound to GTP, focusing on two flexible loops and on the region associated with oncogenic mutations. We perform microsecond-long biased simulations on the wild type and on G12C, G12D, G12 V mutants, describe the resulting free-energy landscapes, and compare our predictions with previous experimental and computational studies. The fast interconversion between open and closed macroscopic states and their similar thermodynamic stabilities are observed. The mutation-induced effects include the alternations of the relative stabilities of different conformational states and the introduction of many microscopic metastable states. Together, our results demonstrate the applicability of the HLDA-based protocol for the conformational sampling of multiple flexible regions in folded proteins.

Mechanistic Understanding from Molecular Dynamics in Pharmaceutical Research 2: Lipid Membrane in Drug Design

We review the use of molecular dynamics (MD) simulation as a drug design tool in the context of the role that the lipid membrane can play in drug action, i.e., the interaction between candidate drug molecules and lipid membranes. In the standard "lock and key" paradigm, only the interaction between the drug and a specific active site of a specific protein is considered; the environment in which the drug acts is, from a biophysical perspective, far more complex than this. The possible mechanisms though which a drug can be designed to tinker with physiological processes are significantly broader than merely fitting to a single active site of a single protein. In this paper, we focus on the role of the lipid membrane, arguably the most important element outside the proteins themselves, as a case study. We discuss work that has been carried out, using MD simulation, concerning the transfection of drugs through membranes that act as biological barriers in the path of the drugs, the behavior of drug molecules within membranes, how their collective behavior can affect the structure and properties of the membrane and, finally, the role lipid membranes, to which the vast majority of drug target proteins are associated, can play in mediating the interaction between drug and target protein. This review paper is the second in a two-part series covering MD simulation as a tool in pharmaceutical research; both are designed as pedagogical review papers aimed at both pharmaceutical scientists interested in exploring how the tool of MD simulation can be applied to their research and computational scientists interested in exploring the possibility of a pharmaceutical context for their research.

Long-Time-Scale Predictions from Short-Trajectory Data: A Benchmark Analysis of the Trp-Cage Miniprotein

Elucidating physical mechanisms with statistical confidence from molecular dynamics simulations can be challenging owing to the many degrees of freedom that contribute to collective motions. To address this issue, we recently introduced a dynamical Galerkin approximation (DGA) [Thiede, E. H. J. Chem. Phys., 150, 2019, 244111], in which chemical kinetic statistics that satisfy equations of dynamical operators are represented by a basis expansion. Here, we reformulate this approach, clarifying (and reducing) the dependence on the choice of lag time. We present a new projection of the reactive current onto collective variables and provide improved estimators for rates and committors. We also present simple procedures for constructing suitable smoothly varying basis functions from arbitrary molecular features. To evaluate estimators and basis sets numerically, we generate and carefully validate a data set of short trajectories for the unfolding and folding of the trp-cage miniprotein, a well-studied system. Our analysis demonstrates a comprehensive strategy for characterizing reaction pathways quantitatively.

Proof of concept for poor inhibitor binding and efficient formation of covalent adducts of KRASG12C and ARS compounds

The use of selective covalent inhibitors with low binding affinity and high reactivity with the target enzyme is a promising way to solve a long-standing problem of the "undruggable" RAS-like proteins. Specifically, compounds of the ARS family that prevent the activation of the GDP-bound G12C mutant of Kirsten RAS (KRAS) are in the focus of recent experimental research. We report the first computational characterization of the entire reaction mechanism of the covalent binding of ARS-853 to the KRASG12C·GDP complex. The application of molecular dynamics, molecular docking and quantum mechanics/molecular mechanics approaches allowed us to model the inhibitor binding to the protein and the chemical reaction of ARS-853 with Cys12 in the enzyme binding site. We estimated a full set of kinetic constants and carried out numerical kinetic analysis of the process. Thus, we were able to compare directly the physicochemical parameters of the reaction obtained in silico and the macroscopic parameters observed in experimental studies. From our computational results, we explain the observed unusual dependence of the rate constant of covalent complex formation, kobs, on the ARS concentration. The latter depends both on the non-covalent binding step with the equilibrium constant, Ki, and on the rate constant of covalent adduct formation, kinact. The calculated ratio kinact/Ki = 213 M-1 s-1 reproduces the corresponding experimental value of 250 ± 30 M-1 s-1 for the interaction of ARS-853 with KRASG12C. Electron density analysis in the reactive region demonstrates that covalent bond formation occurs efficiently according to the Michael addition mechanism, which assumes the activation of the C[double bond, length as m-dash]C bond of ARS-853 by a water molecule and Lys16 in the binding site of KRASG12C. We also refine the kinact and Ki constants of the ARS-107 compound, which shares common features with ARS-853, and show that the decrease in the kinact/Ki ratio in the case of ARS-107 is explained by changes in both Ki and kinact constants.

Insulin Dissociates by Diverse Mechanisms of Coupled Unfolding and Unbinding

The protein hormone insulin exists in various oligomeric forms, and a key step in binding its cellular receptor is dissociation of the dimer. This dissociation process and its corresponding association process have come to serve as paradigms of coupled (un)folding and (un)binding more generally. Despite its fundamental and practical importance, the mechanism of insulin dimer dissociation remains poorly understood. Here, we use molecular dynamics simulations, leveraging recent developments in umbrella sampling, to characterize the energetic and structural features of dissociation in unprecedented detail. We find that the dissociation is inherently multipathway with limiting behaviors corresponding to conformational selection and induced fit, the two prototypical mechanisms of coupled folding and binding. Along one limiting path, the dissociation leads to detachment of the C-terminal segment of the insulin B chain from the protein core, a feature believed to be essential for receptor binding. We simulate IR spectroscopy experiments to aid in interpreting current experiments and identify sites where isotopic labeling can be most effective for distinguishing the contributions of the limiting mechanisms.

SAMPL6 host-guest binding affinities and binding poses from spherical-coordinates-biased simulations

Host-guest binding is a challenging problem in computer simulation. The prediction of binding affinities between hosts and guests is an important part of the statistical assessment of the modeling of proteins and ligands (SAMPL) challenges. In this work, the volume-based variant of well-tempered metadynamics is employed to calculate the binding affinities of the host-guest systems in the SAMPL6 challenge. By biasing the spherical coordinates describing the relative position of the host and the guest, the initial-configuration-induced bias vanishes and all possible binding poses are explored. The agreement between the predictions and the experimental results and the observation of new binding poses indicate that the volume-based technique serves as a nice candidate for the calculation of binding free energies and the search of the binding poses.

Role of the N-Terminal Transmembrane Helix Contacts in the Activation of FGFR3

Fibroblast growth factor receptor 3 (FGFR3) is a member of receptor tyrosine kinases, which is involved in skeletal cell growth, differentiation, and migration. FGFR3 transduces biochemical signals from the extracellular ligand-binding domain to the intracellular kinase domain through the conformational changes of the transmembrane (TM) helix dimer. Here, we apply generalized replica exchange with solute tempering method to wild type (WT) and G380R mutant (G380R) of FGFR3. The dimer interface in G380R is different from WT and the simulation results are in good agreement with the solid-state nuclear magnetic resonance (NMR) spectroscopy. TM helices in G380R are extended more than WT, and thereby, G375 in G380R contacts near the N-termini of the TM helix dimer. Considering that both G380R and G375C show the constitutive activation, the formation of the N-terminal contacts of the TM helices can be generally important for the activation mechanism. © 2019 Wiley Periodicals, Inc.

Quantitative Characterization of Protein-Lipid Interactions by Free Energy Simulation between Binary Bilayers

Using a recently developed binary bilayer system (BBS) consisting of two patches of laterally contacting bilayers, umbrella sampling molecular dynamics (MD) simulations were performed for quantitative characterization of protein-lipid interactions. The BBS is composed of 1,2-dilauroyl-sn-glycero-3-phosphocholine (DLPC) and 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) with an embedded model membrane protein, a gramicidin A (gA) channel. The calculated free energy difference for the transfer of a gA channel from DLPC (hydrophobic thickness ≈ 21.5 Å) to DMPC (hydrophobic thickness ≈ 25.5 Å) bilayers, ΔG(DLPC → DMPC), is -2.2 ± 0.7 kcal/mol. This value appears at odds with the traditional view that the hydrophobic length of the gA channel is ∼22 Å. To understand this discrepancy, we first note that recent MD simulations by different groups have shown that lipid bilayer thickness profiles in the vicinity of a gA channel differ qualitatively from the deformation profile predicted from continuum elastic bilayer models. Our MD simulations at low and high gA:lipid molar ratios and different membrane compositions indicate that the gA channel's effective hydrophobic length is ∼26 Å. Using this effective hydrophobic length, ΔG(DLPC → DMPC) determined here is in excellent agreement with predictions based on continuum elastic models (-3.0 to -2.2 kcal/mol) where the bilayer deformation energy is approximated as a harmonic function of the mismatch between the channel's effective hydrophobic length and the hydrophobic thickness of the bilayer. The free energy profile for gA in the BBS includes a barrier at the interface between the two bilayers which can be attributed to the line tension at the interface between two bilayers with different hydrophobic thicknesses. This observation implies that translation of a peptide between two different regions of a cell membrane (such as between the liquid ordered and disordered phases) may include effects of a barrier at the interface in addition to the relative free energies of the species far from the interface. The BBS allows for direct transfer free energy calculations between bilayers without a need of a reference medium, such as bulk water, and thus provides an efficient simulation protocol for the quantitative characterization of protein-lipid interactions at all-atom resolution.

Replica-Exchange Umbrella Sampling Combined with Gaussian Accelerated Molecular Dynamics for Free-Energy Calculation of Biomolecules

We have developed an enhanced conformational sampling method combining replica-exchange umbrella sampling (REUS) with Gaussian accelerated molecular dynamics (GaMD). REUS enhances the sampling along predefined reaction coordinates, while GaMD accelerates the conformational dynamics by adding a boost potential to the system energy. The method, which we call GaREUS (Gaussian accelerated replica-exchange umbrella sampling), enhances the sampling more efficiently than REUS or GaMD, while the computational resource for GaREUS is the same as that required for REUS. The two-step reweighting procedure using the multistate Bennett acceptance ratio method and the cumulant expansion for the exponential average is applied to the simulation trajectories for obtaining the unbiased free-energy landscapes. We apply GaREUS to the calculations of free-energy landscapes for three different cases: conformational equilibria of N-glycan, folding of chignolin, and conformational change of adenyl kinase. We show that GaREUS speeds up the convergences of free-energy calculations using the same amount of computational resources as REUS. The free-energy landscapes reweighted from the trajectories of GaREUS agree with previously reported ones. GaREUS is applicable to free-energy calculations of various biomolecular dynamics and functions with reasonable computational costs.

U-shaped caveolin-1 conformations are tightly regulated by hydrogen bonds with lipids

The structure and dynamics of a truncated (residues 82-136) caveolin-1 (Cav1) construct having a helix-break-helix motif are explored by both all-atom free energy and molecular dynamics (MD) simulations in an explicit bilayer membrane. Two stable Cav1 conformations with small (LB-Cav1) and large hinge angles (RB-Cav1) between two helices are identified although their relative free energy cannot be reliably estimated due to the sampling issues. RB-Cav1s contain one or two lipids residing between the helices that are hydrogen bonded (h-bonded) to both helices in a multidentate fashion. LB-Cav1s show the helices with mono-dentate lipid h-bond interactions or multidentate interactions limited to a single helix at most. The two conformational states of Cav1 remain their initial state during 2-μs MD simulation, suggesting that there is a significant hidden barrier (other than the insertion depth of Cav1 and its hinge angle) and the Cav1 conformational states are tightly regulated by the h-bonds between Cav1 and lipids along with the associated lipid rearrangement during the course of Cav1 conformational changes. © 2019 Wiley Periodicals, Inc.

Replica-Exchange Methods for Biomolecular Simulations

In this study, a replica-exchange method was developed to overcome conformational sampling difficulties in computer simulations of spin glass or other systems with rugged free-energy landscapes. This method was then applied to the protein-folding problem in combination with molecular dynamics (MD) simulation. Owing to its simplicity and sampling efficiency, the replica-exchange method has been applied to many other biological problems and has been continuously improved. The method has often been combined with other sampling techniques, such as umbrella sampling, free-energy perturbation, metadynamics, and Gaussian accelerated MD (GaMD). In this chapter, we first summarize the original replica-exchange molecular dynamics (REMD) method and discuss how new algorithms related to the original method are implemented to add new features. Heterogeneous and flexible structures of an N-glycan in a solution are simulated as an example of applications by REMD, replica exchange with solute tempering, and GaMD. The sampling efficiency of these methods on the N-glycan system and the convergence of the free-energy changes are compared. REMD simulation protocols and trajectory analysis using the GENESIS software are provided to facilitate the practical use of advanced simulation methods.

A generalized linear response framework for expanded ensemble and replica exchange simulations

Expanded ensemble simulation is a powerful technique for enhancing sampling over a range of thermodynamic parameters. However, although the premise is relatively simple, running successful simulations in practice still presents something of an ad hoc challenge. Three main difficulties exist: (1) the selection of the thermodynamic states, (2) the selection of the sampling weights, and (3) efficient sampling of the expanded parameter space. Here we consider these problems in the context of a pairwise linear response approach to the work fluctuation theorem. The approach offers comprehensive tactics for addressing the three difficulties and can be used as either an alternative or a complement to replica exchange simulations. Importantly, the results are trivially implemented for multi-dimensional parameter spaces and they recover results from the literature aimed at the special cases of simulated/parallel tempering and replica exchange umbrella sampling. Illustrative examples are shown using the NAMD simulation engine.

Quantitative Characterization of Cholesterol Partitioning between Binary Bilayers

We have devised a practical simulation protocol for quantitative characterization of cholesterol (Chol) partitioning between bilayers with different lipid types. The simulation model contains two patches of laterally contacting lipid bilayers, where the host lipids of each bilayer are allowed to self-adjust their packing. For two combinations of bilayers with different lipid types, 1,2-dioleoyl- sn-glycero-3-phosphocholine (DOPC)/1-palmitoyl-2-oleoyl- sn-glycero-3-phosphocholine (POPC) and 1,2-dipalmitoyl- sn-glycero-3-phosphocholine (DPPC), the simulation model has been verified by self-adjusted lipid packing in each bilayer, convergence of Chol partitioning between different Chol initial distributions, and relative diffusion coefficients consistent to those from experiments. The calculated Chol partition coefficient between POPC and DOPC bilayers from the Chol partitioning simulations in the POPC-DPPC and DOPC-DPPC binary bilayer systems shows an excellent agreement with that from available Chol exchange experiments between 1-stearoyl-2-oleoyl- sn-glycero-3-phosphocholine(SOPC)/DOPC vesicles and β-cyclodextrins, which further validates the simulation protocol and illustrates its applicability to any molecular partitioning in the binary bilayer system.

Molecular dynamics simulations of lipid nanodiscs

A lipid nanodisc is a discoidal lipid bilayer stabilized by proteins, peptides, or polymers on its edge. Nanodiscs have two important connections to structural biology. The first is associated with high-density lipoprotein (HDL), a particle with a variety of functionalities including lipid transport. Nascent HDL (nHDL) is a nanodisc stabilized by Apolipoprotein A-I (APOA1). Determining the structure of APOA1 and its mimetic peptides in nanodiscs is crucial to understanding pathologies related to HDL maturation and designing effective therapies. Secondly, nanodiscs offer non-detergent membrane-mimicking environments and greatly facilitate structural studies of membrane proteins. Although seemingly similar, natural and synthetic nanodiscs are different in that nHDL is heterogeneous in size, due to APOA1 elasticity, and gradually matures to become spherical. Synthetic nanodiscs, in contrast, should be homogenous, stable, and size-tunable. This report reviews previous molecular dynamics (MD) simulation studies of nanodiscs and illustrates convergence and accuracy issues using results from new multi-microsecond atomistic MD simulations. These new simulations reveal that APOA1 helices take 10-20 μs to rearrange on the nanodisc, while peptides take 2 μs to migrate from the disc surfaces to the edge. These systems can also become kinetically trapped depending on the initial conditions. For example, APOA1 was trapped in a biologically irrelevant conformation for the duration of a 10 μs trajectory; the peptides were similarly trapped for 5 μs. It therefore remains essential to validate MD simulations of these systems with experiments due to convergence and accuracy issues. This article is part of a Special Issue entitled: Emergence of Complex Behavior in Biomembranes edited by Marjorie Longo.

Free-energy simulations reveal molecular mechanism for functional switch of a DNA helicase

Helicases play key roles in genome maintenance, yet it remains elusive how these enzymes change conformations and how transitions between different conformational states regulate nucleic acid reshaping. Here, we developed a computational technique combining structural bioinformatics approaches and atomic-level free-energy simulations to characterize how the Escherichia coli DNA repair enzyme UvrD changes its conformation at the fork junction to switch its function from unwinding to rezipping DNA. The lowest free-energy path shows that UvrD opens the interface between two domains, allowing the bound ssDNA to escape. The simulation results predict a key metastable 'tilted' state during ssDNA strand switching. By simulating FRET distributions with fluorophores attached to UvrD, we show that the new state is supported quantitatively by single-molecule measurements. The present study deciphers key elements for the 'hyper-helicase' behavior of a mutant and provides an effective framework to characterize directly structure-function relationships in molecular machines.

Dimer Interface of the Human Serotonin Transporter and Effect of the Membrane Composition

The oligomeric state of membrane proteins has recently emerged in many cases as having an effect on their function. However, the intrinsic dynamics of their spatial organization in cells and model systems makes it challenging to characterize. Here we use molecular dynamics (MD) simulations at multiple resolutions to determine the dimer conformation of the human serotonin transporter (hSERT). From self-assembly simulations we predict dimer candidates and subsequently quantify their relative strength. We use umbrella sampling (US) replica exchange MD simulations for which we present extensive analysis of their efficiency and improved sampling compared to regular US MD simulations. The data shows that the most stable hSERT dimer interface is symmetrical and involves transmembrane helix 12 (TM12), similar to the crystal structure of the bacterial homologue LeuT, but with a slightly different orientation. We also describe the supramolecular organization of hSERT from a 250 μs self-assembly simulation. Finally, the effects of the presence of phosphatidylinositol bisphosphate or cholesterol in the membrane model has been quantified for the TM12-TM12 predicted interface. Collectively, the presented data bring new insight to the area of protein and lipid interplay in biological membranes.

Conformational Heterogeneity of the HIV Envelope Glycan Shield

To better understand the conformational properties of the glycan shield covering the surface of the HIV gp120/gp41 envelope (Env) trimer, and how the glycan shield impacts the accessibility of the underlying protein surface, we performed enhanced sampling molecular dynamics (MD) simulations of a model glycosylated HIV Env protein and related systems. Our simulation studies revealed a conformationally heterogeneous glycan shield with a network of glycan-glycan interactions more extensive than those observed to date. We found that partial preorganization of the glycans potentially favors binding by established broadly neutralizing antibodies; omission of several specific glycans could increase the accessibility of other glycans or regions of the protein surface to antibody or CD4 receptor binding; the number of glycans that can potentially interact with known antibodies is larger than that observed in experimental studies; and specific glycan conformations can maximize or minimize interactions with individual antibodies. More broadly, the enhanced sampling MD simulations described here provide a valuable tool to guide the engineering of specific Env glycoforms for HIV vaccine design.

Optimal allosteric stabilization sites using contact stabilization analysis

Proteins can be destabilized by a number of environmental factors such as temperature, pH, and mutation. The ability to subsequently restore function under these conditions by adding small molecule stabilizers, or by introducing disulfide bonds, would be a very powerful tool, but the physical principles that drive this stabilization are not well understood. The first problem lies is in choosing an appropriate binding site or disulfide bond location to best confer stability to the active site and restore function. Here, we present a general framework for predicting which allosteric binding sites correlate with stability in the active site. Using the Karanicolas-Brooks Gō-like model, we examine the dynamics of the enzyme β-glucuronidase using an Umbrella Sampling method to thoroughly sample the conformational landscape. Each intramolecular contact is assigned a score termed a "stabilization factor" that measures its correlation with structural changes in the active site. We have carried out this analysis for three different scaling strengths for the intramolecular contacts, and we examine how the calculated stabilization factors depend on the ensemble of destabilized conformations. We further examine a locally destabilized mutant of β-glucuronidase that has been characterized experimentally, and show that this brings about local changes in the stabilization factors. We find that the proximity to the active site is not sufficient to determine which contacts can confer active site stability. © 2016 Wiley Periodicals, Inc.

Simulating the Activation of Voltage Sensing Domain for a Voltage-Gated Sodium Channel Using Polarizable Force Field

Voltage-gated sodium (Na_V) channels play vital roles in the signal transduction of excitable cells. Upon activation of a Na_V channel, the change of transmembrane voltage triggers conformational change of the voltage sensing domain, which then elicits opening of the pore domain and thus allows an influx of Na⁺ ions. Description of this process with atomistic details is in urgent demand. In this work, we simulated the partial activation process of the voltage sensing domain of a prokaryotic Na_V channel using a polarizable force field. We not only observed the conformational change of the voltage sensing domain from resting to preactive state, but also rigorously estimated the free energy profile along the identified reaction pathway. Comparison with the control simulation using an additive force field indicates that voltage-gating thermodynamics of Na_V channels may be inaccurately described without considering the electrostatic polarization effect.

Interplay of G Protein-Coupled Receptors with the Membrane: Insights from Supra-Atomic Coarse Grain Molecular Dynamics Simulations

G protein-coupled receptors (GPCRs) are central to many fundamental cellular signaling pathways. They transduce signals from the outside to the inside of cells in physiological processes ranging from vision to immune response. It is extremely challenging to look at them individually using conventional experimental techniques. Recently, a pseudo atomistic molecular model has emerged as a valuable tool to access information on GPCRs, more specifically on their interactions with their environment in their native cell membrane and the consequences on their supramolecular organization. This approach uses the Martini coarse grain (CG) model to describe the receptors, lipids, and solvent in molecular dynamics (MD) simulations and in enough detail to allow conserving the chemical specificity of the different molecules. The elimination of unnecessary degrees of freedom has opened up large-scale simulations of the lipid-mediated supramolecular organization of GPCRs. Here, after introducing the Martini CGMD method, we review these studies carried out on various members of the GPCR family, including rhodopsin (visual receptor), opioid receptors, adrenergic receptors, adenosine receptors, dopamine receptor, and sphingosine 1-phosphate receptor. These studies have brought to light an interesting set of novel biophysical principles. The insights range from revealing localized and heterogeneous deformations of the membrane bilayer at the surface of the protein, specific interactions of lipid molecules with individual GPCRs, to the effect of the membrane matrix on global GPCR self-assembly. The review ends with an overview of the lessons learned from the use of the CGMD method, the biophysical-chemical findings on lipid-protein interplay.

Convergence and Sampling in Determining Free Energy Landscapes for Membrane Protein Association

Potential of mean force (PMF) calculations are used to characterize the free energy landscape of protein–lipid and protein–protein association within membranes. Coarse-grained simulations allow binding free energies to be determined with reasonable statistical error. This accuracy relies on defining a good collective variable to describe the binding and unbinding transitions, and upon criteria for assessing the convergence of the simulation toward representative equilibrium sampling. As examples, we calculate protein–lipid binding PMFs for ANT/cardiolipin and Kir2.2/PIP₂, using umbrella sampling on a distance coordinate. These highlight the importance of replica exchange between windows for convergence. The use of two independent sets of simulations, initiated from bound and unbound states, provide strong evidence for simulation convergence. For a model protein–protein interaction within a membrane, center-of-mass distance is shown to be a poor collective variable for describing transmembrane helix–helix dimerization. Instead, we employ an alternative intermolecular distance matrix RMS (D_RMS) coordinate to obtain converged PMFs for the association of the glycophorin transmembrane domain. While the coarse-grained force field gives a reasonable K_d for dimerization, the majority of the bound population is revealed to be in a near-native conformation. Thus, the combination of a refined reaction coordinate with improved sampling reveals previously unnoticed complexities of the dimerization free energy landscape. We propose the use of replica-exchange umbrella sampling starting from different initial conditions as a robust approach for calculation of the binding energies in membrane simulations.

QM/MM free energy simulations: recent progress and challenges

Due to the higher computational cost relative to pure molecular mechanical (MM) simulations, hybrid quantum mechanical/molecular mechanical (QM/MM) free energy simulations particularly require a careful consideration of balancing computational cost and accuracy. Here we review several recent developments in free energy methods most relevant to QM/MM simulations and discuss several topics motivated by these developments using simple but informative examples that involve processes in water. For chemical reactions, we highlight the value of invoking enhanced sampling technique (e.g., replica-exchange) in umbrella sampling calculations and the value of including collective environmental variables (e.g., hydration level) in metadynamics simulations; we also illustrate the sensitivity of string calculations, especially free energy along the path, to various parameters in the computation. Alchemical free energy simulations with a specific thermodynamic cycle are used to probe the effect of including the first solvation shell into the QM region when computing solvation free energies. For cases where high-level QM/MM potential functions are needed, we analyze two different approaches: the QM/MM-MFEP method of Yang and co-workers and perturbative correction to low-level QM/MM free energy results. For the examples analyzed here, both approaches seem productive although care needs to be exercised when analyzing the perturbative corrections.

How Tolerant are Membrane Simulations with Mismatch in Area per Lipid between Leaflets?

Difficulties in estimating the correct number of lipids in each leaflet of complex bilayer membrane simulation systems make it inevitable to introduce a mismatch in lipid packing (i.e., area per lipid) and thus alter the lateral pressure of each leaflet. To investigate potential impacts of such mismatch on simulation results, we performed molecular dynamics simulations of saturated and monounsaturated lipid bilayers with and without gramicidin A or WALP23 at various mismatches by adjusting the number of lipids in the lower leaflet from no mismatch to a 25% reduction compared to that in the upper leaflet. All simulations were stable under the constant pressure barostat, but the mismatch induces asymmetric lipid packing between the leaflets, so that the upper leaflet becomes more ordered, and the lower leaflet becomes less ordered. The mismatch impacts on various bilayer properties are mild up to 5–10% mismatch, and bilayers with fully saturated chains appear to be more prone to these impacts than those with unsaturated tails. The nonvanishing leaflet surface tensions and the free energy derivatives with respect to the bilayer curvature indicate that the bilayer would be energetically unstable in the presence of mismatch. We propose a quantitative criterion for allowable mismatch based on the energetics derived from a continuum elastic model, which grows as a square root of the number of the lipids in the system. On the basis of this criterion, we infer that the area per lipid mismatch up to 5% would be tolerable in various membrane simulations of reasonable all-atom system sizes (40–160 lipids per leaflet).

Coupled folding and binding with 2D Window-Exchange Umbrella Sampling

Intrinsically disordered regions of proteins can gain structure by binding to a partner. This process, of coupled folding and binding (CFaB), is a fundamental part of many important biological processes. Structure-based models have proven themselves capable of revealing fundamental aspects of how CFaB occurs, however, typical methods to enhance the sampling of these transitions, such as replica exchange, do not adequately sample the transition state region of this extremely rare process. Here, we use a variant of Umbrella Sampling to enforce sampling of the transition states of CFaB of HdeA monomers at neutral pH, an extremely rare process that occurs over timescales ranging from seconds to hours. Using high resolution sampling in the transition state region, we cluster states along the principal transition path to obtain a detailed description of coupled binding and folding for the HdeA dimer, revealing new insight into the ensemble of states that are accessible to client recognition. We then demonstrate that exchanges between umbrella sampling windows, as done in previous work, significantly improve relaxation in variables orthogonal to the restraints used. Altogether, these results show that Window-Exchange Umbrella Sampling is a promising approach for systems that exhibit flexible binding, and can reveal transition state ensembles of these systems in high detail. © 2015 Wiley Periodicals, Inc.

Molecular dynamics simulations of biological membranes and membrane proteins using enhanced conformational sampling algorithms

This paper reviews various enhanced conformational sampling methods and explicit/implicit solvent/membrane models, as well as their recent applications to the exploration of the structure and dynamics of membranes and membrane proteins. Molecular dynamics simulations have become an essential tool to investigate biological problems, and their success relies on proper molecular models together with efficient conformational sampling methods. The implicit representation of solvent/membrane environments is reasonable approximation to the explicit all-atom models, considering the balance between computational cost and simulation accuracy. Implicit models can be easily combined with replica-exchange molecular dynamics methods to explore a wider conformational space of a protein. Other molecular models and enhanced conformational sampling methods are also briefly discussed. As application examples, we introduce recent simulation studies of glycophorin A, phospholamban, amyloid precursor protein, and mixed lipid bilayers and discuss the accuracy and efficiency of each simulation model and method. This article is part of a Special Issue entitled: Membrane Proteins edited by J.C. Gumbart and Sergei Noskov.

Simulations of outer membrane channels and their permeability

Channels in the outer membrane of Gram-negative bacteria provide essential pathways for the controlled and unidirectional transport of ions, nutrients and metabolites into the cell. At the same time the outer membrane serves as a physical barrier for the penetration of noxious substances such as antibiotics into the bacteria. Most antibiotics have to pass through these membrane channels to either reach cytoplasmic bound targets or to further cross the hydrophobic inner membrane. Considering the pharmaceutical significance of antibiotics, understanding the functional role and mechanism of these channels is of fundamental importance in developing strategies to design new drugs with enhanced permeation abilities. Due to the biological complexity of membrane channels and experimental limitations, computer simulations have proven to be a powerful tool to investigate the structure, dynamics and interactions of membrane channels. Considerable progress has been made in computer simulations of membrane channels during the last decade. The goal of this review is to provide an overview of the computational techniques and their roles in modeling the transport across outer membrane channels. A special emphasis is put on all-atom molecular dynamics simulations employed to better understand the transport of molecules. Moreover, recent molecular simulations of ion, substrate and antibiotics translocation through membrane pores are briefly summarized. This article is part of a Special Issue entitled: Membrane Proteins edited by J.C. Gumbart and Sergei Noskov.

Combine umbrella sampling with integrated tempering method for efficient and accurate calculation of free energy changes of complex energy surface

Umbrella sampling is an efficient method for the calculation of free energy changes of a system along well-defined reaction coordinates. However, when there exist multiple parallel channels along the reaction coordinate or hidden barriers in directions perpendicular to the reaction coordinate, it is difficult for conventional umbrella sampling to reach convergent sampling within limited simulation time. Here, we propose an approach to combine umbrella sampling with the integrated tempering sampling method. The umbrella sampling method is applied to chemically more relevant degrees of freedom that possess significant barriers. The integrated tempering sampling method is used to facilitate the sampling of other degrees of freedom which may possess statistically non-negligible barriers. The combined method is applied to two model systems, butane and ACE-NME molecules, and shows significantly improved sampling efficiencies as compared to standalone conventional umbrella sampling or integrated tempering sampling approaches. Further analyses suggest that the enhanced performance of the new method come from the complemented advantages of umbrella sampling with a well-defined reaction coordinate and integrated tempering sampling in orthogonal space. Therefore, the combined approach could be useful in the simulation of biomolecular processes, which often involves sampling of complex rugged energy landscapes.

Lipid-linked oligosaccharides in membranes sample conformations that facilitate binding to oligosaccharyltransferase

Lipid-linked oligosaccharides (LLOs) are the substrates of oligosaccharyltransferase (OST), the enzyme that catalyzes the en bloc transfer of the oligosaccharide onto the acceptor asparagine of nascent proteins during the process of N-glycosylation. To explore LLOs' preferred location, orientation, structure, and dynamics in membrane bilayers of three different lipid types (dilauroylphosphatidylcholine, dimyristoylphosphatidylcholine, and dioleoylphosphatidylcholine), we have modeled and simulated both eukaryotic (Glc3-Man9-GlcNAc2-PP-Dolichol) and bacterial (Glc1-GalNAc5-Bac1-PP-Undecaprenol) LLOs, which are composed of an isoprenoid moiety and an oligosaccharide, linked by pyrophosphate. The simulations show no strong impact of different bilayer hydrophobic thicknesses on the overall orientation, structure, and dynamics of the isoprenoid moiety and the oligosaccharide. The pyrophosphate group stays in the bilayer head group region. The isoprenoid moiety shows high flexibility inside the bilayer hydrophobic core, suggesting its potential role as a tentacle to search for OST. The oligosaccharide conformation and dynamics are similar to those in solution, but there are preferred interactions between the oligosaccharide and the bilayer interface, which leads to LLO sugar orientations parallel to the bilayer surface. Molecular docking of the bacterial LLO to a bacterial OST suggests that such orientations can enhance binding of LLOs to OST.

GENESIS: a hybrid-parallel and multi-scale molecular dynamics simulator with enhanced sampling algorithms for biomolecular and cellular simulations

GENESIS (Generalized-Ensemble Simulation System) is a new software package for molecular dynamics (MD) simulations of macromolecules. It has two MD simulators, called ATDYN and SPDYN. ATDYN is parallelized based on an atomic decomposition algorithm for the simulations of all-atom force-field models as well as coarse-grained Go-like models. SPDYN is highly parallelized based on a domain decomposition scheme, allowing large-scale MD simulations on supercomputers. Hybrid schemes combining OpenMP and MPI are used in both simulators to target modern multicore computer architectures. Key advantages of GENESIS are (1) the highly parallel performance of SPDYN for very large biological systems consisting of more than one million atoms and (2) the availability of various REMD algorithms (T-REMD, REUS, multi-dimensional REMD for both all-atom and Go-like models under the NVT, NPT, NPAT, and NPγT ensembles). The former is achieved by a combination of the midpoint cell method and the efficient three-dimensional Fast Fourier Transform algorithm, where the domain decomposition space is shared in real-space and reciprocal-space calculations. Other features in SPDYN, such as avoiding concurrent memory access, reducing communication times, and usage of parallel input/output files, also contribute to the performance. We show the REMD simulation results of a mixed (POPC/DMPC) lipid bilayer as a real application using GENESIS. GENESIS is released as free software under the GPLv2 licence and can be easily modified for the development of new algorithms and molecular models. WIREs Comput Mol Sci 2015, 5:310–323. doi: 10.1002/wcms.1220

Cholesterol modulates the dimer interface of the β₂-adrenergic receptor via cholesterol occupancy sites

The β2-adrenergic receptor is an important member of the G-protein-coupled receptor (GPCR) superfamily, whose stability and function are modulated by membrane cholesterol. The recent high-resolution crystal structure of the β2-adrenergic receptor revealed the presence of possible cholesterol-binding sites in the receptor. However, the functional relevance of cholesterol binding to the receptor remains unexplored. We used MARTINI coarse-grained molecular-dynamics simulations to explore dimerization of the β2-adrenergic receptor in lipid bilayers containing cholesterol. A novel (to our knowledge) aspect of our results is that receptor dimerization is modulated by membrane cholesterol. We show that cholesterol binds to transmembrane helix IV, and cholesterol occupancy at this site restricts its involvement at the dimer interface. With increasing cholesterol concentration, an increased presence of transmembrane helices I and II, but a reduced presence of transmembrane helix IV, is observed at the dimer interface. To our knowledge, this study is one of the first to explore the correlation between cholesterol occupancy and GPCR organization. Our results indicate that dimer plasticity is relevant not just as an organizational principle but also as a subtle regulatory principle for GPCR function. We believe these results constitute an important step toward designing better drugs for GPCR dimer targets.

Applications of rare event dynamics on the free energy calculations for membrane protein systems

Techniques of rare event dynamics were reviewed including the string methods, which will be implemented with the biochemical simulation packages. The existing methods were applied to study biological systems with relevance to drug design and drug metabolism. The rare event dynamics simulations were performed to understand the kinetic and thermodynamic free energy information on the drug binding sites in the M2 proton channel, and also the free energy of insertion and association of membrane proteins and membrane active peptides. Results give a theoretical framework to interpret and reconcile existing and often conflicting opinions.

Model for coupled insertion and folding of membrane-spanning proteins

Current understanding of the forces directing the folding of integral membrane proteins is very limited compared to the detailed picture available for water-soluble proteins. While mechanistic studies of the folding process in vitro have been conducted for only a small number of membrane proteins, the available evidence indicates that their folding process is thermodynamically driven like that of soluble proteins. In vivo, however, the majority of integral membrane proteins are installed in membranes by dedicated machinery, suggesting that the cellular systems may act to facilitate and regulate the spontaneous physical process of folding. Both the in vitro folding process and the in vivo pathway must navigate an energy landscape dominated by the energetically favorable burial of hydrophobic segments in the membrane interior and the opposition to folding due to the need for passage of polar segments across the membrane. This manuscript describes a simple, exactly solvable model which incorporates these essential features of membrane protein folding. The model is used to compare the folding time under conditions which depict both the in vitro and in vivo pathways. It is proposed that the cellular complexes responsible for insertion of membrane proteins act by lowering the energy barrier for passage of polar regions through the membrane, thereby allowing the chain to more rapidly achieve the folded state.

Theory of Adaptive Optimization for Umbrella Sampling

We present a theory of adaptive optimization for umbrella sampling. With the analytical bias force constant obtained from the constrained thermodynamic length along the reaction coordinate, the windows are distributed to optimize the overlap between neighbors. Combining with the replica exchange method, we propose a method of adaptive window exchange umbrella sampling. The efficiency gain in sampling by the present method originates from the optimal window distribution, in which windows are concentrated to the region where the free energy is steep, as well as consequently improved random walk.

Computational Recipe for Efficient Description of Large-Scale Conformational Changes in Biomolecular Systems

Characterizing large-scale structural transitions in biomolecular systems poses major technical challenges to both experimental and computational approaches. On the computational side, efficient sampling of the configuration space along the transition pathway remains the most daunting challenge. Recognizing this issue, we introduce a knowledge-based computational approach toward describing large-scale conformational transitions using (i) nonequilibrium, driven simulations combined with work measurements and (ii) free energy calculations using empirically optimized biasing protocols. The first part is based on designing mechanistically relevant, system-specific reaction coordinates whose usefulness and applicability in inducing the transition of interest are examined using knowledge-based, qualitative assessments along with nonequilirbrium work measurements which provide an empirical framework for optimizing the biasing protocol. The second part employs the optimized biasing protocol resulting from the first part to initiate free energy calculations and characterize the transition quantitatively. Using a biasing protocol fine-tuned to a particular transition not only improves the accuracy of the resulting free energies but also speeds up the convergence. The efficiency of the sampling will be assessed by employing dimensionality reduction techniques to help detect possible flaws and provide potential improvements in the design of the biasing protocol. Structural transition of a membrane transporter will be used as an example to illustrate the workings of the proposed approach.

Multidimensional umbrella sampling and replica-exchange molecular dynamics simulations for structure prediction of transmembrane helix dimers

Structural information of a transmembrane (TM) helix dimer is useful in understanding molecular mechanisms of important biological phenomena such as signal transduction across the cell membrane. Here, we describe an umbrella sampling (US) scheme for predicting the structure of a TM helix dimer in implicit membrane using the interhelical crossing angle and the TM-TM relative rotation angles as the reaction coordinates. This scheme conducts an efficient conformational search on TM-TM contact interfaces, and its robustness is tested by predicting the structures of glycophorin A (GpA) and receptor tyrosine kinase EphA1 (EphA1) TM dimers. The nuclear magnetic resonance (NMR) structures of both proteins correspond to the global free-energy minimum states in their free-energy landscapes. In addition, using the landscape of GpA as a reference, we also examine the protocols of temperature replica-exchange molecular dynamics (REMD) simulations for structure prediction of TM helix dimers in implicit membrane. A wide temperature range in REMD simulations, for example, 250-1000 K, is required to efficiently obtain a free-energy landscape consistent with the US simulations. The interhelical crossing angle and the TM-TM relative rotation angles can be used as reaction coordinates in multidimensional US and be good measures for conformational sampling of REMD simulations.